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Preface 


The  purpose  of  the  research  presented  in  this  paper  was  to  compare 
several  control  variates  for  queueing  network  simulation.  During  the 
literature  search  and  review  I  noted  the  most  promising  internal  and 
external  control  variates.  I  also  included  a  new  external  control 
variate  which  was  calculated  using  Bell  Lab's  Queueing  Network  Analyzer, 
a  network  decomposition  algorithm  employing  two-moment  approximations  of 
the  stochastic  processes  in  the  network. 

The  experiment  was  designed  to  be  general  enough  to  apply  to  many 
queueing  networks.  The  variance  ratios  of  each  of  the  control  variates 
provides  a  measure  of  their  efficiency  to  reduce  the  variance  of  the 
parameter  being  estimated.  The  coverage  probabilities  of  the  confidence 
intervals  formed  about  the  controlled  estimates  of  the  mean  responses 
provide  a  measure  of  the  accuracy  of  the  veriance  reduction  technique. 
That  is,  a  control  variate  that  results  in  a  confidence  interval  with 
poor  coverage  of  the  true  parameter  is  not  very  useful  even  if  the 
variance  reduction  is  significant. 

I  wish  to  thank  my  thesis  advisor,  Major  Joseph  R.  Litko,  Ph.D., 
and  Major  Kenneth  W.  Bauer,  Jr.,  Ph.D.  for  their  guidance  in  this 
research  effort.  I  also  wish  to  thank  Major  Anthony  P.  Sharon  for 
providing  me  with  a  copy  of  his  thesis  which  was  the  foundation  and 
inspiration  of  this  work.  Finally,  I  wish  to  thank  my  wife  Jamie  for 
her  understanding  and  support  during  the  many  days  and  nights  I  was  tied 
to  my  computer. 

John  J.  Tomick 
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Abstract 


The  purpose  of  this^-s-tudy  mss  to  compare  several  control  variates 
for  queueing  network  simulation.  The  author's  goal  was  to  provide  the 
simulation  community  with  some  guidance  for  selecting  control  variates 
that  will  lead  to  significant  reductions  in  the  variance  of  the 
estimated  responses  that  do  not  introduce  significant  bias. 

Both  internal  and  external  control  variates  were  examined.  The 
measures  for  comparing  them  were  the  variance  ratios  obtained  for  each 
control  variate  against  each  of  the  two  response  variables  and  the 
coverage  of  their  respective  9SZ  confidence  intervals. 

The  two  response  variables  selected  for  this  research  were  the 
average  sojourn  time  in  the  network  and  the  probability  that  the  number 
in  the  fourth  queue  exceeds  twice  the  mean  number  in  queue  at  steady- 
state.  The  internal  controls  included  standardized  routing  controls  and 
standardized  work  variables.  The  external  controls  included  the  average 
sojourn  time  in  the  network  and  the  average  number  in  queue  at  the 
fourth  node.  The  external  controls  were  further  classified  into  two 
groups:  analytic  Jackson  controls  and  analytic  approximations.  K  r  ) 
The  analytic  Jackson  control  variates  were  found  by  decomposing 
the  network  and  using  the  M/M/1  formulas.  The  analytic  approximations 
were  found  using  C/C/1  formulas,  specifically  those  employed  in  Bell 
Laboratories'  Queueing  Network  Analyzer.  The  observed  values  of  the 
external  control  variates  were  found  by  using  the  following  parameters 
observed  during  a  run  of  the  simulation  model:  the  external  arrival 
rate,  the  mean  service  times  at  each  node,  the  squared  coefficient  of 


x 


variation  of  the  service  tines  at  each  node,  and  the  probabilistic 
routing  matrix.  The  "known"  neans  of  the  external  control  variates  were 
found  using  the  values  of  the  above  parameters  that  were  input  to  the 
simulation  model. 

The  network  of  queues  studied  was  an  open  network  with  Poisson 
external  arrivals  and  exponential  servers  at  the  first  three  nodes.  At 
the  fourth  node  the  service  times  were  generated  from  the  maximum 
entropy  and  the  hyperexponential  distributions  characterized  by  two 
moments. 

In  general,  the  external  control  variates  achieved  smaller 
variance  ratios  than  the  internal  control  variates;  however,  the 
coverages  of  the  confidence  intervals  about  the  controlled  responses 
were  worse.  The  range  of  the  average  variance  ratios  for  the  internal 
control  variates  was  0.886  to  0.949  with  coverages  from  0.898  to  0.929 
for  the  95X  confidence  Intervals.  The  range  of  the  average  variance 
ratios  for  the  external  control  variates  was  0.494  to  0.774  with 
coverages  from  0.576  to  0.775  for  the  95X  confidence  intervals. 
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A  COMPARISON  OF  CONTROL  VARIATES 
FOR  QUEUEING  NETWORK  SIMULATION 


I.  Introduction 

The  purpose  of  this  paper  Is  to  report  the  results  of  the  author's 
research  comparing  the  effectiveness  and  the  bias  of  several  control 
variates  for  queueing  network  simulation.  For  this  research,  the  author 
selected  some  of  the  most  promising  control  variates  found  in  the 
literature.  Also,  the  author  included  a  new  external  control  technique 
that  makes  use  of  a  software  package  which  approximates  the  performance 
measures  of  a  queueing  network. 

The  inspiration  and  foundation  for  this  research  came  from  the 
work  done  by  Sharon  (1986).  He  examined  the  effectiveness  of  Jackson 
networks  as  control  variates  for  queueing  network  simulation.  A  more 
complete  discussion  of  Sharon's  work  can  be  found  in  the  review  of  the 
literature  presented  in  Chapter  II. 

The  author's  research  makes  a  significant  contribution  to  the 
experiential  knowledge  of  control  variates  in  two  areas.  First,  most  of 
the  literature  on  experimental  results  of  control  variates  report  only 
the  efficiency  of  the  technique  in  terms  of  variance  ratios,  percent 
variance  reduction,  etc.  This  research  effort  also  examined  the  bias 
introduced  to  the  estimates  through  the  use  of  control  variates  in  terms 
of  the  estimated  coverage  of  the  95Z  confidence  intervals  about  the 
controlled  responses.  Second,  this  research  is  the  first  attempt  to 
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compare  both  internal  and  external  control  variates.  The  ultimate  goal 
of  this  research  was  to  provide  the  simulation  community  with  some 
guidance  for  selecting  control  variates  for  use  in  queueing  network 
simulation. 

Following  the  literature  review,  the  author  discusses  the 
methodology  of  his  approach  to  the  research  in  Chapter  III.  Then,  the 
results  of  the  research  are  presented  in  Chapter  IV.  Finally,  the 
author's  conclusions  and  recommendations  are  given  in  Chapter  V. 

The  rest  of  this  chapter  is  devoted  to  Introducing  the  reader  to 
some  of  the  basic  terms,  concepts,  and  notation  used  in  queueing  theory 
and  in  simulation  modeling,  respectively.  The  reader  who  is  familiar 
with  these  areas  may  skip  to  Chapter  II. 

Queueing  Theory 

Basic  queueing  theory  involves  customers  arriving  to  a  service 
center.  The  customers  can  represent  people  waiting  in  line  for  a  bank 
teller,  cars  waiting  in  line  at  a  toll  booth,  and  so  forth.  Because 
queueing  theory  has  such  broad  applicability,  one  will  usually  see  the 
the  more  generic  term  entity  used  instead  of  customer. 

A  Single  Queue.  To  describe  a  queueing  system  we  need  a 
description  of  the  arrival  process,  a  description  of  the  service 
mechanism,  and  a  queue  discipline. 

The  arrival  process  can  be  deterministic;  that  is,  the  time 
between  arrivals  of  entitles  is  a  constant.  It  can  also  be  a  random 
process  described  by  a  probability  density  function.  Traditionally,  the 
arrival  rate  (the  mean  number  of  entitles  arriving  to  a  node  per  unit 
time)  is  denoted  by  the  Greek  letter  lambda  (X).  The  mean  interarrival 
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tiae,  or  aean  time  between  arrivals,  is  slaply  the  reciprocal  of  the 
arrival  rate. 

Similarly,  the  service  aechanisa  uy  be  deterainistlc  or 
probabilistic.  Furtheraore,  the  service  center  aay  have  aore  than  one 
identical  server.  Traditionally,  the  service  rate  (the  aean  nuaber  of 
entities  that  can  be  serviced  per  unit  tine)  Is  denoted  by  the  Creek 
letter  au  (p).  Likewise,  the  aean  service  tine,  or  aean  tiae  between 
service  coapletlons,  is  the  reciprocal  of  the  service  rate. 

The  traffic  intensity  at  a  node  is  defined  to  be  the  ratio  of  the 
arrival  rate  to  the  product  of  the  service  rate  and  the  nuaber  of 
identical  servers,  which  is  denoted  by  the  Creek  letter  rho  (/»). 
Matheaatically,  p  *  X/bji  ,  where  a  denotes  the  nuaber  of  identical 
servers.  If  the  traffic  intensity  is  greater  than  one,  then  entities 
are  arriving  to  the  queue  faster  than  they  can  be  serviced.  Therefore, 
the  queue  grows  infinitely  long  (i.e.  the  queue  is  unstable)  unless  the 
calling  population  is  finite  or  the  queue  is  capacitated  with  blocking 
or  balking.  Uhen  a  queue  has  a  finite  capacity  with  blocking,  then  the 
servers  that  feed  the  queue  stop  servicing  entitles  until  there  is  rooa 
for  thea  in  the  following  queue.  On  the  other  hand,  a  queue  filled  to 
capacity  that  allows  balking  causes  entitles  arriving  to  the  queue  to  be 
routed  to  another  queue  or  to  leave  the  systea  entirely.  The  research 
presented  in  this  paper  exaained  a  network  of  queues  with  infinite 
capacities  at  equllibriua  (or  steady  state),  which  iaplies  a  traffic 
Intensity  less  than  unity. 

The  traffic  intensity  is  also  referred  to  as  the  utilization 
factor  for  the  service  center,  since  it  represents  the  fraction  of  the 
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server's  capacity  that  la  being  utilized  on  the  average  by  arriving 
entitles. 


Finally,  the  queue  discipline  indicates  how  to  select  the  next 
entity  waiting  for  service.  The  nost  couon  queue  discipline  is  first 
in.  first  out  (FIFO),  which  is  also  referred  to  as  first  come,  first 
served  (FCFS).  However,  one  night  also  specify  last  in,  first  out 
(LIFO)  as  the  queue  discipline,  or  set  up  soae  sort  of  priority 
selection  based  upon  the  type  of  entity. 

A  Hetwork  of  Queues.  Networks  of  queues  are  useful  nodels  for 
describing  many  real-world  systems,  such  as  computer  time-sharing 
processes,  communication  systems,  transportation  systems,  and  assembly¬ 
line  operations. 

A  queueing  network  is  composed  of  nodes  representing  a  system's 
service  mechanisms  and  directed  arcs  between  the  nodes  representing  the 
flow  of  entities  through  the  system.  The  individual  nodes  are 
conventionally  labeled  by  the  distribution  of  the  interarrival  tines, 
the  distribution  of  the  service  times,  and  the  number  of  identical 
servers.  For  example,  an  M/M/1  queue  has  exponentially  distributed 
interarrival  times  and  service  times  with  one  server.  The  "H"  stands 
for  Markovian  (memoryless),  and  the  exponential  distribution  is  the  only 
continuous  distribution  with  the  Markovian  property.  In  general,  a 
queue  with  many  Independent  input  processes  that  are  not  necessarily 
Markovian,  a  non-Markovian  service-time  distribution,  and  many  servers 
is  denoted  by  Cl/C/m. 

There  are  two  broad  classifications  of  queueing  networks:  open, 
and  closed.  An  open  queueing  network  allows  for  external  arrivals  to 


the  nodes  (referred  to  as  exogenous  arrivals).  In  general,  entities  nay 
enter  or  leave  the  network  Iron  any  node.  Alternatively,  a  closed 
network  has  a  fixed  nuaber  of  entities  cycling  through  it;  no  additional 
entities  arrive  to  the  network,  and  none  of  the  entities  in  the  network 
ever  leave. 

Jackson  networks.  A  Jackson  network  is  an  analytically  tractable 
queueing  network  nodel.  Although,  it  is  a  fairly  restrictive  nodel,  it 
has  found  many  uses  and  is  a  good  approxiaation  to  aany  real-world 
systeas.  In  a  Jackson  network,  all  external  arrivals  are  Independent 
Poisson  processes  (that  is,  their  interarrival  tiaes  are  independent  and 
exponentially  distributed),  all  servers  have  exponentially  distributed 
service  tines,  the  queue  discipline  is  FIFO,  the  queue  capacity  at  each 
node  is  infinite,  the  routing  between  nodes  can  be  probabilistic  but  not 
conditional,  and  the  tine  to  travel  between  two  nodes  is  zero. 

Burke's  theoren  states  that  "the  steady-state  output  of  a  stable 
M/H/n  queue  with  input  paraneter  A  and  service-tine  parameter  \i  for  each 
of  the  a  channels  is  in  fact  a  Poisson  process  at  the  sane  rate  A" 
(Kleinrock,  1975:149).  Therefore,  given  a  network  of  M/M/a  queues,  the 
input  process  to  any  node  (i)  in  the  network  is  a  Poisson  process  with 
paraneter  A,,  which  is  a  aixture  of  the  output  processes  of  nodes 
feeding  into  it  plus  the  external  arrival  process.  And,  as  Jackson 
discovered  "each  node  ...  in  the  network  behaves  as  if  it  were  an 
independent  N/N/a  systen  with  a  Poisson  input  rate  A,”  (Kleinrock, 
1975:150).  This  allows  a  network  to  be  decoaposed  into  independent 
H/H/a  queues  whose  perforaance  measures  can  be  solved  for  analytically. 
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If  the  queueing  network  can  be  represented  as  a  Jackson  network, 
then  the  performance  neasures  (such  as  the  average  waiting  tise  at  a 
particular  queue,  the  utilization  of  a  particular  server,  or  the  average 
tine  an  entity  spends  in  the  systen)  can  be  solved  for  analytically.  If 
we  drop  the  subscript  (i)  for  readability,  then  the  steady-state  results 
for  the  N/N/n  queue  at  any  node  can  be  solved  for  using  the  following 
equations: 


P  »  A/np 


P. 


P. 


v  (A/p)*  <A/|1>-  1 

2,  - + - 

ni  nl  (1  -  p) 

<A/p)* 

- - —  P.,  if  0  1  n  i  i 

nl 

(A/p)* 

■-  ...  P..  if  nil 

nln 


P.(A/p)'p 
*  *  «l(l  -  /»)’ 

W,  =  L,/A 
L  =  L,  +  (A/p) 
W  .  W,  +  (1/p) 


where 

A  s  arrival  rate  of  entities  to  the  node 
n  =  number  of  identical  servers 
p  *  service  rate  of  servers  at  the  node 
p  «  traffic  intensity  or  server  utilization 
P«  *  probability  that  there  are  0  entities  in  the  system 
P„  >  probability  that  there  are  n  entities  in  the  systen 
L  >  expected  number  of  entities  in  the  systen 
L,  •  expected  queue  length  (excludes  entitles  in  service) 
H  a  expected  time  in  system 
V,  a  expected  waiting  time  (excludes  service  tine) 


(1.1) 

(1.2) 


(1.3) 


(1.4) 

(1.5) 

(1.6) 
(1.7) 
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However,  there  are  aany  systems  for  which  the  assumptions  of  the 
Jackson  network  are  grossly  violated.  This  leads  to  a  queueing  model 
that  is  analytically  Intractable.  The  performance  measures  of  such 
networks  can  be  solved  for  by  methods  which  use  approximations,  or  by 
simulating  the  model.  The  reader  is  referred  to  Chapter  II  for  the 
author's  review  of  analytical  approaches  to  solving  queueing  networks 
using  approximations. 

Simulation  Modeling 

Simulation  is  an  experimental  technique  for  analyzing  complex 
systems  which  are  analytically  Intractable  and  usually  Involve 
stochastic  processes  (indexed  collections  of  random  variables).  As 
such,  the  output  from  a  simulation  experiment  is  a  random  variable;  and 
therefore,  the  measured  response  is  only  an  estimate  of  the  true 
parameter  of  interest.  The  variance  of  the  estimate  is  a  measure  of  its 
precision;  and,  in  most  cases,  variance  reduction  techniques  provide  a 
means  of  obtaining  more  precise  estimates  with  minimal  cost  in  terms  of 
computer  resources. 

Variance  Reduction  Techniques.  In  a  recent  survey  of  variance 
reduction  techniques  (VRTs),  Wilson  classified  "...  all  VRTs  into  two 
major  categories — correlation  methods  and  Importance  methods"  (Wilson, 
1984:280).  In  the  paper  he  discusses  three  correlation  methods  (common 
random  numbers,  antithetic  variates,  and  control  variates)  and  four 
Importance  methods  (importance  sampling,  conditional  Honte  Carlo, 
stratified  sampling,  and  systematic  sampling). 

The  basic  distinction  between  the  two  categories  is  the  underlying 
principle  of  the  techniques.  Correlation  methods  increase  the 
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efficiency  of  a  simulation  by  making  use  of  the  linear  correlations 
among  simulation  responses,  and  importance  methods  use  prior  knowledge 
of  the  input  domain  to  achieve  a  variance  reduction.  For  a  more 
complete  discussion  of  these  techniques  the  reader  is  referred  to 
Kleljnen  (1974),  Law  and  Kelton  (1982),  and  Wilson  (1984). 

Sharon's  research  and  the  research  presented  in  this  paper  use 
control  variates  to  achieve  a  variance  reduction  in  the  estimates  of  the 
performance  measures  of  interest.  A  control  variate  must  have  a  known 
expectation  and  be  correlated  with  the  response. 

Methods  for  Obtaining  Control  Variates.  Law  and  Kelton  describe 
three  general  methods  for  obtaining  control  variates.  The  first  method 
uses  the  correlation  between  the  input  random  varlable(s)  and  the  output 
random  varlable(s).  Since  these  kinds  of  control  variates  must  be 
generated  during  the  simulation  they  are  called  internal  or  concomi  tant 
control  variates. 

A  second  method  involves  the  simulation  of  a  similar  system  that 
is  analytically  tractable.  Using  common  random  numbers,  the 
corresponding  output  of  the  second  simulation  becomes  the  control 
variate  for  the  system  under  study.  These  kinds  of  control  variates  are 
called  external  control  variates.  Note  that  this  method  assumes  that 
there  is  a  significant  correlation  between  the  results  of  the  two 
systems  through  the  use  of  common  random  numbers. 

The  third  method  makes  use  of  the  situations  when  there  are 


several  unbiased  estimators  of  the  mean  of  the  performance  measure  of 
interest.  In  such  cases  a  new  controlled  estimator  can  be  formed  as  a 
convex  combination  of  the  existing  estimators  (Law  and  Kelton,  1982:358- 
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II .  Literature  Review 

The  following  discussion  Is  a  review  of  the  literature  that  Is 
applicable  to  the  author's  research  Into  the  following  two  areas:  (1) 
control  variates  for  queueing  network  simulation,  and  (2)  two-moment 
approximations  to  performance  measures  of  the  GI/G/m  queue.  The 
discussion  is  presented  in  a  topical  format  as  outlined  below. 

A.  Control  Variates 

1.  Justification  for  New  Research 

2.  Theory  of  Control  Variates 

3.  Results  of  Previous  Research 

B.  Approximations  to  Point  Processes 

1.  Need  for  Approximations 

2.  Two-Moment  Approximations 

3.  The  Queueing  Network  Analyzer 

Control  Variates 

In  the  first  chapter,  the  author  introduced  variance  reduction 
techniques  which  are  used  to  provide  more  precise  estimates  of  the 
response  from  a  simulation  experiment.  This  research  compared  and 
contrasted  several  control  variates  for  their  efficiency  in  reducing  the 
variance  of  the  estimated  responses  and  for  the  amount  of  bias  that  may 
have  been  introduced  to  the  estimates. 

Justification  for  New  Research.  As  mentioned  previously,  there 
are  two  major  categories  of  variance  reduction  techniques— correlation 
methods  and  importance  methods.  Of  all  the  techniques,  "...  the  method 
of  control  variates  is  one  of  the  most  promising"  (Lavenberg  and  Welch, 
1981:322). 
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Correlation  Methods.  The  correlation  methods  include  coaaon 
randoa  numbers,  antithetic  variates,  and  control  variates.  A  brief 
discussion  of  the  uses  and  drawbacks  of  these  correlation  aethods  are 
presented  below. 

Coaaon  Randoa  Huabers.  The  aethod  of  coaaon  randoa 
nuabers  can  be  used  when  comparing  two  or  aore  alternative  systea 
designs,  or  when  designing  an  experiaent  for  a  response  surface  model. 
This  aethod  assuaes  the  existence  of  a  positive  correlation  between  the 
randoa  number  streams  driving  the  simulation  and  the  measured 
response! s) . 

However,  in  complex  simulation  models,  especially  queueing  network 
models,  the  correlation  tends  to  be  very  weak.  In  soae  experiments  it 
aay  even  be  negative  resulting  in  a  variance  increase.  Furthermore,  the 
synchronization  of  coaaon  randoa  number  streams  aay  be  difficult  if  not 
impossible  for  soae  simulation  experiaents. 

Antithetic  Variates.  The  aethod  of  antithetic 
variates  is  applicable  to  the  simulation  of  a  single  system.  This 
aethod  tries  to  Induce  a  negative  correlation  between  pairs  of  runs  of 
the  simulation  aodel  to  achieve  a  variance  reduction.  According  to  Law 
and  Kelton,  the  method  "...  dates  back  at  least  to  1956  with  the  paper 
of  Haanersley  and  Norton  in  the  context  of  Monte  Carlo  simulation"  (Law 
and  Kelton,  1982:354).  More  recently,  Schruben  and  Hargolln  (1978) 
demonstrated  the  effectiveness  of  antithetic  variates  in  conjunction 
with  a  2*  factorial  experimental  design. 


But,  this  aethod  suffers  froa  the  sane  drawbacks  as  does  the 
aethod  of  coaaon  randoa  nuabers — the  correlation  may  be  weak  or  opposite 


in  sign  to  that  desired,  and  the  synchronization  of  random  number 
streams  may  be  a  problem.  Also,  there  are  some  other  assumptions  that 
must  be  met  to  use  Schruben  and  Margolin's  assignment  rule.  As  with 
common  random  numbers,  the  method  of  antithetic  variates  seems  to  be 
more  successful  when  used  with  Monte  Carlo  simulation  experiments. 

Control  Variates.  The  method  of  control  variates  is 
applicable  to  any  simulation  experiment  Involving  stochastic  processes 
with  known  means.  A  practitioner  may  simultaneously  collect 
observations  for  several  internal  controls  and  select  those  having  a 
significant  correlation  with  the  response  variable(s).  Furthermore,  the 
practitioner  can  calculate  the  magnitude  of  the  variance  reduction 
without  further  simulation  runs.  If  the  practitioner  had  used  common 
random  numbers  or  antithetic  variates  and  wanted  to  know  by  how  much  the 
variances  of  the  estimates  were  reduced,  it  would  require  additional 
simulation  runs  using  independent  random  number  streams. 

However,  control  variates  have  some  drawbacks  as  well.  The 
traditional  type  of  external  control  variates  require  additional 
simulation  runs  and,  the  resulting  controlled  estimator  may  be  biased. 

Importance  Methods.  The  importance  methods  include  several 
sampling  techniques,  which  have  not  found  much  popularity  among 
simulation  practitioners.  Furthermore,  Pritsker  concluded  that  the 
Importance  methods  require  further  refinement  before  they  can  be  applied 
to  complex  simulation  experiments  (Pritsker,  1986:749). 

Theory  of  Control  Variates.  The  fundamental  idea  behind  the 
method  of  control  variates  is  to  select  a  random  variable  with  a  known 
expectation  that  is  highly  correlated  with  the  response  variable. 
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Let  Y 


Univariate  Simulation  Response  with  a  Single  Control, 
be  an  unbiased  estimator  of  the  parameter  of  interest  •  ;  that  is,  the 
expectation  of  Y,  denoted  E(Y),  equals  •  .  Let  C  be  another  random 
variable  with  known  expectation  pe  that  is  highly  correlated  with  Y. 
Then,  for  any  constant  b  (known  as  the  control  coefficient),  the 
controlled  estiaator  Y(b),  given  by  Eq  (2.1),  is  unbiased  for  ft. 

Y(b)  *  Y  -  b(C  -  |ie)  (2.1) 

The  variance  of  Y(b)  is  given  by 

VarlY(b) 1  *  Var(Y)  +  b*Var(C)  -  2bCov(Y,C)  (2.2) 

and  a  variance  reduction  will  be  realized  if 

2bCov(Y,C)  >  b*Var(C)  (2.3) 

That  is,  if  Eq  (2.3)  is  satisfied,  then  the  controlled  estiaator  will 
have  a  smaller  variance  than  the  uncontrolled  estimator.  With  a  little 
calculus  it  is  easy  to  show  froa  Eq  (2.2)  that  Y(b)  has  ainlaua  variance 
when  b  is  set  equal  to  the  optiaal  control  coefficient,  £,  which  is 
given  by 


0  -  Cov(Y,C)/Var(C)  (2.4) 

Substituting  Eq  (2.4)  into  Eq  (2.1)  leads  to  the  optimal 
controlled  estimator  Y(0),  which  is  given  by 

Y(0)  .  Y  -  CCov(Y,C)/Var(C) ]• (C  -  p.)  (2.5) 


2.4 


Anderson  (1984)  provides  a  proof  that  the  variance  of  the  controlled 
estimator  which  is  given  by 

Var(Y(0) ]  *  (1  -  Pre')‘Va r(Y)  (2.6) 

where  p Te*  is  the  square  of  the  correlation  coefficient  between  the 
response  variable  Y  and  and  the  control  variate  C.  Because  the 
correlation  coefficient  in  Eq  (2.6)  is  squared,  the  sign  of  the 
correlation  does  not  matter;  only  the  size  of  the  correlation  is 
important.  As  I e I  ■»  1,  the  correlation  between  Y  and  C  becomes  more 
significant,  and  the  size  of  the  variance  reduction  increases. 

Let  8,  the  parameter  of  interest,  be  denoted  by  |iT.  Then,  we  know 
that  the  average  of  the  uncontrolled  observations  Y,  is  an  unbiased 
point  estimator  of  pT.  Furthermore,  the  average  of  the  controlled 
observations  Y,(0)  is  also  an  unbiased  estimator  of  uT.  This  is 
represented  mathematically  as  follows: 

_  K 

Y(0)  »  (l/K)  l  (2.7) 

I  •  I 

where  K  is  the  sample  size  and 

Y,(0)  *  Y,  -  0(C,  -  pe)  (2.8) 

In  practice,  Cov(Y,C)  and  Var(C)  are  unknown;  and  therefore,  6  is 
unknown  and  must  be  estimated.  Following  Bauer  (1987),  the  intuitive 
approach  to  estimating  0  is  to  replace  the  right-hand  side  of  Eq  (2.4) 
with  the  appropriate  sample  statistics,  which  yields  the  least-squares 
solution.  Under  the  assumption  of  joint  normality  between  Y  and  C,  the 
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least  squares  solution  is  also  the  maximum  likelihood  solution  (Bauer, 
1987:6).  Then,  (3  can  be  estimated  by 

.  <  _  * 

e  =>  I<Y,  -  Y)(C,  -  C)/  £<C,  -  C)*  (2.9) 

I  •  1  I  •  I 

where 

Y  =  1  Y,/R  (2.10) 

t  •  I 

and 

K 

C  =  l  C./K  (2.11) 

I  •  I 

Furthermore,  a  point  estimate  of  yT  is  given  by 

« 

Y(0)  »  XY,(0)/K  (2.12) 

<  •  t 

or 

Y(0)  «  Y  -  0(C  -  Me)  (2.13) 

And,  the  variance  of  the  point  estimator  is  given  by 

Var[Y<en  -  Var[Y(0)]/R  (2.14) 

where 

Vart Y(8) 1  =  (1  -  p,rc)-Var(Y)  (2.15) 

Bauer  (1987)  provides  the  derivation  of  the  interval  estimate 
under  the  assumption  that  Y  and  C  are  jointly  normal  random  variates. 
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The  resulting  100(l-a)X  confidence  interval  is  given  by  the  following 
equation 


Y<0)  ±  t« .*(l-o/2)  •  (Var[Y(|8)  ]*s, ,} 1  '*  (2.16) 

where 

R  K 

s,,  »  X(C,  -  M,)*/K  £(C,  -  C)*  (2.17) 

1  ■  I  I  •  t 

and  t,.,(l-«/2)  is  the  100(l-a/2)  percentile  of  Student's  t-distribution 
with  (K-2)  degrees  of  freedoa. 

Since  0  Bust  be  estimated,  we  expect  to  achieve  a  smaller  variance 
reduction  than  that  which  could  have  been  obtained  had  we  known  the 
optimal  control  coefficient.  Lavenberg,  Moeller  and  Welch  (1982) 
quantified  this  loss  by  what  is  known  as  the  loss  factor  (LF).  It  is 
defined  as  "the  ratio  of  the  variance  of  the  estimator  of  |iT  when  the 
optimal  control  coefficient  is  not  known  to  the  variance  of  the 
estimator  when  the  coefficient  is  known"  (Bauer,  1987:9).  Bauer 
provides  the  derivation  of  the  loss  factor,  which  reduces  to 

LF  -  (I-2)/(K-Q-2)  (2.18) 

where 

Q  ■  the  number  of  controls  (for  the  univariate  case  Q=l) 

This  "loss  factor  acts  as  a  multiplier  to  the  minimum  variance  ratio 
(MVR)"  (Bauer,  1987:10,14)  given  by 

MVR  >  7ar[Y(0)  J/Var(Y)  (2.19) 
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The  MVS  represents  the  possible  variance  reduction  when  the  optical 
control  coefficient  is  known.  Multiplying  Eq  (2.18)  and  Eq  (2. 19) 
together  leads  to  the  variance  ratio  (VR).  The  VR  represents  the 
possible  variance  reduction  when  0  is  not  known. 

VR  -  LF-MVR  (2.20) 

Univariate  Simulation  Response  with  Multiple  Controls. 
Kleijnen  (1974)  addresses  the  extensions  to  multiple  control  variates. 
Also,  Bauer  (1987)  provides  a  summary  of  "the  development  presented  by 
Lavenberg  and  Welch  (1981)  for  simulation  output  analysis  based  on 
independent  replications,  batch  means,  and  regenerative  analysis" 

(Bauer,  1987:11). 

Let  Y  be  the  univariate  response  with  variance  ar* ,  C  be  the  (QX1) 

aw 

vector  of  controls,  tfCT  be  the  (QX1)  vector  of  covariances  between  Y  and 

mt 

C,  and  E«  be  the  (QXQ)  covariance  matrix  of  the  controls.  Then, 
rewriting  Eq  (2.13)  with  multiple  controls  leads  to 


Y((9)  *  Y  -  £t(C  -  a.)  (2.21) 

Mr  mt  rv  Mr 

A  _ 

where  0,  C,  and  ji«  are  (QX1)  vectors.  The  vector  of  optimal  control 
coefficients,  is  then  given  by 

0  *  E,“«r,T  (2.22) 

MT  MT  MT 

Since  the  covariance  matrices  are  usually  unknown,  0  can  be  estimated  by 
substituting  the  sample  analogs  of  E,  and  tfer  into  Eq  (2.22).  This 
leads  to  the  following  equation: 
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0  9  S*  Sgi 


(2.23) 


where  S«"‘  is  the  inverse  of  the  (QXQ)  saeple  covariance  matrix  of  the 
controls,  and  SCr  is  the  (QX1)  vector  of  sample  covariances  between  the 
univariate  response  and  the  vector  of  controls. 

0nder  the  assumption  that  Y  and  £  have  the  Joint  multivariate 
normal  distribution 


(2.2A) 


Y(|J)  is  unbiased  for  (iv,  and  an  exact  100(l-a)Z  confidence  interval  is 
given  by 


Y<8)  +  t«.,.,(l-a/2)D-ST.e  (2.25) 

where 

D*  *  r‘  ♦  (*-l)“(C  -  m,)tSc*‘(C  -  u.)  (2.26) 

M  IV  M  fSt  Jw 

s».,*  -  (I-Q-D  ’d-lXS,*  -  S,vfS.*'Stf)  (2.27) 

m  mr  iv 

t«.t.,(l-a/2)  is  the  100(l-a/2)  percentile  of  Student's  t -distribution 
with  (t-Q-l)  degrees  of  freedom,  and  STs  is  the  sample  variance  of  Y 
(Bauer  and  others,  1988:3-4).  Experimental  results  have  shown  that  the 
assumption  of  joint  multivariate  normality  is  robust  (Bauer,  1988). 

Multiple  Simulation  Responses  with  Multiple  Controls. 

Bauer,  Venkatraman  and  Wilson  (1987)  provide  an  outline  of  the 
theoretical  formulas  for  the  case  when  there  are  P  response  variables 
and  Q  control  variables.  In  terms  of  the  notation,  the  univariate 
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response  Y  becomes  a  (PX1)  vector  of  response  variables,  0  becomes  a 
( PXQ)  matrix  of  control  coefficients,  and  the  scaler  sample  standard 
deviation  S ,  becomes  the  sample  covariance  matrix  of  the  response 
vector.  Under  the  assumption  that  Y  and  C  have  the  joint  multivariate 

M  M 

normal  distribution 


(2.28) 


Y(0)  is  an  unbiased  estimator  of  ja T ,  and  an  exact  100(l-a)X  confidence 
ellipsoid  for  ji,  is  given  by 


*  ww)(w-Q)"D-r(i-«;p,i-p-Q)  (2.29) 

where 

D*  -  1“  ♦  U-1)*‘(C  -  M«)TSe*‘(C  -  |i.)  (2.30) 

**  Jw  A*  «  « 

s».e*  «  (I-Q-1)*'(*-1)(ST  -  STBSe-'Set)  (2.31) 

A#  M  MM  AT 

and  F(l-a;m,,nt)  is  the  100(l-a)  percentile  of  the  F-distribution  with 
m,  and  m,  degrees  of  freedom  (Bauer  and  others,  1987:335). 

The  advantage  of  the  above  approach  over  selecting  separate 
controls  for  each  response  is  the  capability  to  form  a  joint  confidence 
region  for  the  response  vector,  rather  than  being  limited  to  univariate 
confidence  intervals. 

Results  of  Previous  Research.  The  following  discussion  of  the 
experimental  results  found  in  the  literature  is  presented  in 
chronological  order. 
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Review  of  Cheng  (1978).  Cheng  (1978)  provided 
Interpretations  of  statistically  well-known  formulas  used  in  ordinary 
regression  analysis  to  control  variates  in  siaulation  under  the 
assumption  of  nornallty. 

The  reader  should  note  that  there  is  an  error  in  Equation  (3)  of 
Cheng's  article.  The  correct  equation,  which  appeared  in  Cheng  and 
Feast  (1980),  reads  as  follows: 


«  »  X  -  e’(X  -  p)  (2.32) 

Review  of  Chens  and  Feast  (1980).  Cheng  and  Feast  (1980) 
made  the  statement  that  "practically  all  control  variables  suggested  in 
the  literature  are  of  the  fora  where  their  wean  \i  is  known,  but  ...  [the 
covariance  aatrix]  is  not"  (Cheng  and  Feast,  1980:51).  However,  aore 
recent  literature  follows  their  suggestion  of  using  standardized  subs 
for  control  variables  rather  than  using  saaple  means  or  straight  subs. 
These  standardized  controls  are  of  the  fora 


C 


»  £  T|/(H) ‘ ' * 

i  ■ 


(2.33) 


If  the  Y,'s  are  a-independent  (X,  and  X, , ,  are  independent  if  j  ^  i), 
with  zero  Bean  and  E(X,X, T)  finite,  then  in  the  liait  £  is  normally 
distributed  with  mean  zero  and  covariance  matrix 


B(X,X, T  ♦  X.X,*  +  X,X,T  ♦  ...  +  X,X.T  +  X.X, T ) 


(2.34) 


(Cheng  and  Feast,  1980:52) 
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Review  of  La venture.  Moeller,  and  Welch  (1982).  Lavenberg, 
Moeller,  and  Welch  experimented  with  three  internal  control  variates  in 
closed  queueing  networks  with  more  than  one  class  of  entitles,  where  a 
generic  class  is  denoted  by  the  letter  (d).  Briefly,  the  control 
variates  used  were  work  variables  (the  sum  of  service  times  for  type  d 
entities  per  type  d  event),  flow  variables  (the  fraction  of  type  d 
events  at  node  1),  and  service-time  variables  (the  sample  service  times 
for  type  d  entities  at  node  i). 

They  reported  achieving  the  largest  variance  reductions  using  work 
variables  and  limited  their  research  accordingly.  They  also  reported 
that  as  the  server  utilization  increased  so  did  the  size  of  the  variance 
reduction  in  waiting  time.  The  range  of  the  estimated  variance  ratios 
using  six  work  variables  for  controls  was  from  0.30  to  0.85.  These 
results  translate  to  minimum  variance  ratios  in  the  range  of  0.16  to 
0.77,  which  are  obtained  by  dividing  the  variance  ratios  by  the 
theoretical  loss  factor. 

Review  of  Wilson  and  Pritsker  (1984).  Wilson  and  Pritsker 
experimented  with  poststratif led  sampling  and  standardized  work 
variables  as  variance  reduction  techniques  adapted  to  the  estimation 
methods  of  replication  analysis  (independent  replications  of  a 
simulation)  and  regenerative  analysis  (independent  cycles  within  a 
simulation).  The  standardized  work  variables  are  given  by  Eq  (2.36), 
and  poststratif led  sampling  refers  to  an  importance  technique  which 
groups  the  responses  into  strata  according  to  a  stratification  variate. 
Wilson  and  Pritsker  reported  the  following  reductions  in  the  variance  of 
the  point-estimators  and  in  the  width  of  the  90X  confidence  intervals 
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that  can  be  obtained  with  each  procedure  for  several  closed  and  nixed 
machine-repair  systems: 


Table  2.1.  Experimental  Results  of  Wilson  and  Pritsker  (1984) 


Variance  Reduction  Variance  Confidence- Interval 

Technique  Reduction  Width  Reduction 


Poststratification  10  -  4QX  1  -  20X 

Work  Variables  20  -  90X  10  -  70X 


Review  of  Sharon  (1986).  Sharon  "investigated  two  types  of 
Jackson  control  variates,  external  and  analytic,  for  estimating  the 
utilization  factors  and  waiting  times  in  three  different  queueing 
networks"  (Sharon,  1986 : 35 ) .  He  obtained  control  variates  for  each 
network  using  three  different  service  time  distributions  (exponential, 
Weibull,  and  uniform)  and  two  different  traffic  intensities  (0.5  and 
0.9),  which  resulted  in  eighteen  experiments. 

External  Jackson  Control  Variates.  To  obtain  external 
control  variates  he  used  a  Jackson  network  approximation  to  each  of  the 
more  general  networks  he  studied.  He  substituted  the  exponential 
distribution  for  the  service  time  distributions  with  the  identical  mean 
service  times  used  in  the  original  model.  Of  course,  this  required  a 
second  simulation  for  each  of  the  three  networks.  Then,  the  output  of 
these  second  simulations  were  contrasted  with  the  results  derived 


analytically  to  regress  out  some  of  the  variance  of  the  estimates 
obtained  from  the  simulation  of  the  original  model. 


For  estimating  the  means  of  the  two  response  variables  at  any 
given  node,  the  controlled  observations  were  of  the  fora  given  by  Eq 
(2.8).  For  example,  in  terns  of  estiaating  the  aean  waiting  tine  at  any 
node  in  the  network 

«,(b)  ■  W,  -  b(C,  -  |ie)  (2.35) 

where 

H((b)  =  a  controlled  observation  of  the  average  waiting  tine 
W,  =  an  uncontrolled  observation  of  the  average  waiting  time 
b  =  the  control  coefficient 

C,  *  an  observation  of  the  average  waiting  tine  from  a  second 
siaulation  model  of  the  approximating  Jackson  Network 
Me  s  the  aean  waiting  tine,  which  is  the  analytic  solution  of 
an  H/H/m  queue  given  by  Eq  (l.i),  using  the  mean  arrival 
rates,  aean  service  rates,  and  probabilistic  routing 
structure  that  was  input  to  the  simulation  model 

The  controlled  observations  of  the  utilization  factors  were  defined  in  a 
similar  manner. 

Analytic  Jackson  Control  Variates.  What  Sharon  terms 
an  analytic  Jackson  control  variate  is  "an  amalgam  of  the  internal  and 
external  approaches"  (Sharon,  1986:18).  Instead  of  using  the  input 
random  variables  directly  as  control  variates,  Sharon  substituted  the 
known  means  and  the  observed  averages  of  the  input  random  variables  into 
the  steady-state  equations  for  an  M/H/m  queue  to  derive  the  desired 
control  variates.  The  input  random  variables  included  the  arrival 
rates,  the  service  rates,  and  the  probabilistic  routing  structure. 

Looking  back  to  Eq  (2.35),  only  the  definition  of  C,  changes.  It 
is  now  the  analytic  solution  of  an  M/M/m  queue  using  the  average  arrival 
rates,  average  service  rates,  and  the  probabilistic  routing  observed 
from  the  original  siaulation. 
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Unfortunately,  as  Sharon  indicated,  there  are  two  drawbacks  to 
using  his  analytic  Jackson  control  variates.  First,  he  was  able  to 
obtain  sample  values  of  the  arrival  rates,  the  service  rates,  and 
probabilistic  routing  at  the  high  setting  of  the  traffic  intensity  (0.9) 
such  that  the  derived  traffic  Intensity  using  these  observed  values  was 
greater  than  unity  (Sharon,  1986:32). 

The  second  drawback  of  the  analytic  Jackson  control  variate  is 
that  the  use  of  “the  observed  mean  arrival  rates  and  service  rates  in 
the  Jackson  model  equations  will  result  in  a  biased  control  variate” 
(Sharon,  1986:32-33).  The  severity  of  this  bias  was  left  to  future 
research. 

Results.  Sharon's  results  show  favorable  variance 
reductions  in  the  estimates  of  the  server  utilization  factors  in  the 
range  of  68  to  99  percent.  The  higher  variance  reductions  were  achieved 
at  the  lower  setting  of  the  traffic  intensity  (0.5).  However,  the 
Jackson  analytic  control  variates  for  the  waiting  times  produced  little 
or  no  variance  reduction,  and  in  some  cases,  they  produced  variance 
increases  in  the  estimates  of  the  waiting  times  (Sharon,  1986:73-75). 

Review  of  Bauer.  Venkatraman  and  Wilson  (1987).  Bauer , 
Venkatraman  and  Wilson  report  a  new  control  variate  estimator  which 
makes  use  of  the  cases  when  the  covariance  matrix  of  the  controls  is 
known.  For  the  experiment,  they  selected  the  standardized  work 
variables  given  by  Wilson  and  Pritsker  (1984)  and  the  standardized 
routing  variables  defined  by  Bauer  (1987).  Both  types  of  internal 
control  variates  mentioned  above  have  known  means  and  known  covariance 
matrices. 
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Standardized  Work  Variables.  If  we  assume  that  the 
service  time  process  at  node  j  Is  given  by  the  Independent  and 
Identically  distributed  (IID)  sequence  (D,(j):  l  £  1),  j  *  l,  ...,  g, 
and  we  define  f ,  to  be  the  number  of  service  times  completed  at  node  j 
in  the  period  [0,t],  then  a  standardized  work  variable  for  node  j  is 

*  i 

«,  *  l  [0, ( J )  -  pj/tf,  (2.36) 

ft  •  ft 


where  «,  is  the  frequency  with  which  an  entity  visits  node  j  and  f  is 
the  sum  of  the  f,*s  (Bauer  and  others,  1987:337). 

Standardized  Routing  Variables.  Define  N,(t)  to  be 
the  number  of  entitles  exiting  from  node  J  in  the  time  interval  [0,t]. 
Define  plk  to  be  the  probability  that  an  entity  exiting  from  node  j  will 
go  to  node  k,  and  define  the  indicator  variable  Illk  as  follows: 


I 


<  j  i 


if  the  ith  entity  leaving  node  j  goes  to  node  k 
otherwise 


Then  a  standardized  routing  variable  for  node  j  is  given  by 


R, 


H  I 

*  Z 

>  •  i 


Pi » )/[H, ( t) •  ( 1  -  PiJPi.J*'* 


(2.37) 


Results.  The  new  control  variate  estimator  yielded  a 
confidence-region  that  is  somewhat  larger  than  the  confidence  region 
obtained  using  the  usual  controlled  confidence-region  estimator.  But, 
it  also  demonstrated  more  reliable  coverage  properties.  However,  the 
100(1-«)X  confidence  ellipsoid  for  pT  is  only  approximate.  Further 
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research  is  belag  conducted  to  develop  a  more  refined  estimator  of  the 
covariance  matrix  of  the  controlled  response  vector  and  an  improved 
confidence-region  estimator  for  |if. 

Approximations  to  Point  Processes 

The  author  investigated  the  efficiency  of  applying  two-aoment 
approx i nations  to  control  variates  in  a  manner  sinilar  to  Sharon's 
analytical  controls.  In  effect,  this  technique  uses  approximations  of 
performance  measures  of  the  exact  network  as  external  control  variates. 
In  the  case  of  Sharon's  analytical  controls,  his  technique  uses  exact 
performance  measures  of  an  approximating  network  as  external  control 
variates. 

Heed  for  Approximations.  In  terms  of  queueing  networks  the 
arrival  and  departure  processes  are  point  processes.  If  the  arrival 
processes  are  renewal  processes,  then  the  congestion  measures  of  the 
individual  queues  and  the  entire  network  can  be  solved  for  analytically. 
In  most  cases  the  departure  process  of  a  queue  is  not  a  renewal  process. 
And  since  the  departure  process  of  one  queue  becomes  an  arrival  process 
to  the  next  queue  in  the  network,  then  the  congestion  measures  cannot  be 
solved  for  using  exact  analytic  methods. 

Whitt  (1982)  investigated  simple  approximations  for  stochastic 
point  processes.  He  considered  point  processes  on  the  positive  real 
line  for  which  the  "total  number  of  points  is  infinite  but  the  number  of 
points  in  any  bounded  interval  is  finite"  (Whitt,  1982:129).  Followin* 
Whitt  (1982),  let 

S.  =  the  position  of  the  nth  point  from  the  origin,  n  >  0, 
and  S,  =  0 
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X.  *  S.  -  S..,,  n  >  1  (the  tiae  Interval  between  successive 
points) 

H(t)  *  aaxin  >0:  S,  i  t),  t  £  0  (the  counting  process 
recording  the  nuaber  of  points  in  the  interval  <0,t] 

Then,  "the  stochastic  processes  (Sa>,  (X.),  and  (N(t))  are  three 
different  representations  of  the  sane  point  process'*  (Whitt,  1902:129). 

Let  the  points  in  the  point  process  represent  the  occurrence  of  a 
certain  event,  then  H(t)  represents  the  nuaber  of  tiaes  the  event 
occurred  In  the  tiae  Interval  (0,t).  If  the  tine  intervals  between  each 
occurrence  of  the  event  are  Independent  and  identically  distributed, 
then  (M(t)>  is  called  a  renewal  counting  process.  The  aost  coaaon 
renewal  counting  process  is  the  Poisson  process.  Furtheraore,  (X.)  is 
called  a  renewal  process,  and  if  (N(t)l  is  a  Poisson  process  with  rate 
A,  then  the  X.'s  are  distributed  exponentially  with  aean  1/A. 

However,  if  the  point  process  (M(t)>  is  not  a  renewal  counting 
process,  then  the  interval  sequence  (X.)  is  not  a  renewal  process;  and 
therefore,  the  aodel  becoaes  analytically  intractable.  Under  such 
clrcuastances  aost  practitioners  will  siaulate  the  aodel.  However, 
another  approach  is  to  approxiaate  the  point  process  by  a  renewal 
process  and  then  solve  the  aodel  analytically. 

Two-Hoaent  Approx i nations.  Whitt  (1982)  refers  to  several  authors 

who 


...  suggest  approxiaating  all  the  flows  (point  processes)  in 
a  network  of  queues  by  renewal  processes  characterized  by 
two  paraaeters.  It  was  discovered  that  one  paraaeter 
(representing  the  rate  of  the  process)  is  usually  not  good 
enough,  but  two  paraaeters  (representing  the  rate  and  the 
variability)  often  are  sufficient  (Whitt,  1982:126). 
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Review  of  Whitt  (1982).  Whitt  describes  ways  to  approximate 
a  single  point  process  by  a  renewal  process  in  two  steps:  "first, 
properties  of  the  point  process  are  used  to  specify  a  few  moments  of  the 
Interval  between  renewals;  then  a  convenient  distribution  is  fit  to 
these  moments'*  (Whitt,  1982:125).  There  are  other  appropriate 
parameters,  but  the  parameters  that  he  has  focused  on  are  "the  moments 
of  the  renewal  interval  in  the  approximating  renewal  process”  (Whitt, 
1982:126).  In  the  paper  he  outlines  two  methods  for  specifying  the 
first  few  moments  of  the  renewal  interval — the  stationary-interval 
method  and  the  asymptotic  method.  Briefly, 


the  stationary-interval  method  equates  the  moments  of  the 
renewal  interval  with  the  moments  of  the  stationary  interval 
in  the  point  process  to  be  approximated.  The  asymptotic 
method,  in  an  attempt  to  account  for  the  dependence  among 
successive  intervals,  determines  the  moments  of  the  renewal 
interval  by  matching  the  asymptotic  behavior  of  the  moments 
of  the  sums  of  successive  intervals  (Whitt,  1982:125). 


Review  of  Whitt  (1984).  Later,  in  1984,  Whitt  published 
methods  for  approximating  the  departure  process  of  a  single-server 
queue.  This  result  is  significant  because  the  departure  process  of  one 
queue  becomes  the  arrival  process  to  the  next  queue  in  a  network. 

Whitt  discusses  how  to  use  the  two  methods  above  for  approximating 
the  departure  process.  An  interesting  note  is  that 


the  asymptotic-method  approximation  for  the  departure 
process  is  just  the  arrival  process,  provided  that  the 
arrival  process  is  in  the  class  of  approximating  processes, 
e.g.  a  renewal  process.  Otherwise,  the  approximating 
process  for  [the  departure  process]  obtained  by  the 
asymptotic  method  is  the  same  as  the  approximating  process 
for  [the  arrival  process]  (Whitt,  1984:502). 
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Finally,  Whitt  and  his  associates 


indicated  three  ways  the  appr ox i nations  sight  be  improved: 
(1)  using  the  third  moment,  (2)  using  the  lag-1  correlation 
. ..,  and  (3)  developing  a  hybrid  procedure  ...  However 
[they]  examined  the  last  two  methods  and  did  not  find  an 
improvement  (Whitt,  1984:516). 


Review  of  Albln  and  Rai  (1986).  Albln  and  Rai  studied  two 
queues  in  series  "to  identify  a  renewal  process  to  approximate  the 
departure  process  of  a  EGI,/M/1  queue"  (Albin  and  Rai,  1986:132).  The 
arrival  process  to  the  first  queue  was  a  superposition  of  independent 
stationary  renewal  processes.  Each  queue  had  a  single  server  with 
exponentially  distributed  service  times,  an  infinite  queue  capacity,  and 
FIFO  queue  discipline  (Albin  and  Rai,  1986:130).  A  hybrid  method  of  two 
basic  methods  (the  Poisson  and  the  asymptotic)  led  to  an  average 
absolute  error  in  hybrid  approximations  of  the  expected  number  in  the 
second  queue  of  6Z  compared  to  the  22-41%  error  in  the  basic  methods 
(Albin  and  Rai,  1986:131).  The  Queueing  Network  Analyzer,  discussed 
later  in  this  chapter,  uses  the  stationary-interval  method  to  identify 
renewal -process  approximations  for  departure  processes.  Albin  and  Rai's 
hybrid  method  applies  to  queues  with  exponentially  distributed  service 
times. 

The  squared  coefficient  of  variation  for  approximating  the  renewal 
departure  intervals  using  the  hybrid  method  is  given  by 


c«*  =  wc.’  +  <1  -  w)c,* 


(2.38) 


where  the  squared  coefficient  of  variation  for  approximating  the  renewal 
departure  intervals  using  the  asymptotic  method  (c.1)  and  the  squared 
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coefficient  of  variation  for  approximating  the  renewal  departure 
intervals  using  the  Poisson  method  (c,1)  are  given  by 


c ,  *  *  E(A,/A)c,* 

Ml 


(2.39) 


c,' 


(2.40) 


(A  is  the  total  arrival  rate)  and  the  weighting  coefficient  u  = 

“(n*. Pi, Pi)  is  given  by 

«<*T, />.,/>*)  =  [A  +  Bn*(l  -  /»,)*  +  C/(  1  -  #>,)*]“  (2.41) 


The  symbols  in  Eq  (2.41)  are  defined  as  follows: 


The  coefficients  (A,B,C)  equal 

(1.0,  2.0,  0.05)  for  calculating  the  expected 
number  of  entities  at  the  second  node, 

(1.0,  1.0,  0.05)  for  calculating  the  standard 
deviation  of  the  number  of  entities  at  the 
second  node,  and 

(1.7,  2.3,  0.04)  for  calculating  the  probability 
of  an  entity  being  delayed  at  the  second  node; 


/>,  and  /»,  are  the  traffic  intensities  at  the  respective  nodes;  and 


V 


l  (A, /A) * 


is  the  effective  number  of  component  arrival 
processes. 


The  weighting  function  is  the  result  of  an  experimental  design 
involving  27  combinations  of  pt,  pt,  and  t|*,  where  the  constants  A,  B, 
and  C  were  identified  using  multiple  linear  regression.  "Different 
weighting  functions  are  needed  for  different  congestion  measures  because 
of  the  properties  of  the  basic  methods"  (Albin  and  Kai,  1986:138). 
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This  hybrid  method  works  well  when  the  squared  coefficient  of  variation 
of  each  of  the  interval  renewal  processes  input  to  the  first  queue  are 
in  the  range  (0,9]  (Albin  and  Kai,  1986:  138). 

Review  of  Kimura  (1986).  Kimura  reported  on  a  two-aoaent 
approximation  that  yields  better  results  than  those  achieved  by  the 
Queueing  Network  Analyzer  (which  is  described  later  in  this  chapter). 
Let  EW(GI/G/m)  denote  the  steady-state  Bean  waiting  time  (until 
beginning  service)  in  the  Gl/C/m  queue.  The  approximation  formula  he 
gives  in  his  paper  is 


EW(CI/C/m)  *  (c.%c.* 


)/ 


EW(D/M/m) 


EW(M/D/m) 


2(c.*+c.*-l) 

EW(M/M/m) 


(2.42) 


where  c.*  and  c.*  are  the  coefficients  of  variation  of  the  interarrival 
times  and  service  times  respectively;  /»  is  the  traffic  intensity;  and 
EW(D/M/m) ,  EW(H/0/m)  and  EU(N/H/m)  are  the  steady-state  mean  waiting 
times  in  the  respective  queueing  systems.  Kimura 's  approximation 
formula  given  by  Eq  (2.42)  "is  a  weighted  harmonic  mean  of  the  expected 
waiting  times  for  the  0/N/m,  M/D/m  and  N/H/m  queues  and  it  is  exact  for 
these  queueing  systems'*  (Kimura,  1986:751). 

Note  that  EU(M/M/m)  is  equivalent  to  W,  given  by  Eq  (1.5)  in  the 
first  chapter.  Because  of  some  numerical  difficulties  in  solving  for 
the  exact  analytic  solutions  of  BU(D/M/m)  and  EW(M/D/m),  Kimura  suggests 
using  the  following  approximations: 

EU(D/M/a)  %  ( l/2)[l-4C(m,p) ]*exp[-2(l-p)/3p]*EW(N/N/m)  (2.43) 
EH( N/0/m)  *  (l/2)[l+C(m,p)]-EW(M/H/m)  (2.44) 

where 
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C(m,p)  -  ( 1-p) (a-1 )£ (4+5*) ' ' *  -  2]/(16a?) 


(2.45) 


Xiaura  found  "that  these  approx last ions  are  fairly  accurate  unless  /*  is 
close  to  zero"  (Xiaura,  1986:761). 

Substituting  the  approxiaations  given  by  Eqs  (2.43)  through  (2.45) 
into  Eq  (2.42)  yields  a  siapler  approxiaatlon  that  is  aore  tractable. 
This  siapler  approxiaatlon  is  given  by 


EW(GI/G/a)  »  (l/2)(c.*  +  c.*)k.EW(M/M/a)  (2.46) 


where  the  correction  factor  k  s  k(GI/G/a)  is  defined  by 


and 


k(GI/C/a) 


k(D/M/a) 


1  -  c,» 
k(N/D/a) 


+  c. 


(2.47) 


k(0/N/a)  -  aax{(l-4C(a,/>)  ]*exp[-2( l-/»)/3^J,  10“>  (2.48) 

k(H/D/a)  -  1  +  C(m,/»)  (2.49) 


The  aaxiaua  in  Eq  (2.48)  Is  used  to  avoid  dividing  by  zero  or 
aeanlngless  approxiaations  with  negative  values  (Xiaura,  1986:761). 

The  Queueing  Network  Analyzer.  The  Queueing  Hetwork  Analyzer 
(QHA)  is  a  commercially  available  "software  package  developed  at  Bell 
Laboratories  to  calculate  approximate  congestion  measures  for  a  network 
of  queues"  (Whitt,  1983a: 2779).  The  first  version  of  QUA  operates  under 
the  following  assumptions: 


Assumption  1.  The  network  is  open  rather  than  closed. 
Customers  come  from  outside,  receive  service  at  one  or  more 
nodes,  and  eventually  leave  the  systea. 

Assumption  2.  There  are  no  capacity  contraints.  There  is 
no  Halt  on  the  number  of  customers  that  can  be  in  the 
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entire  network  and  each  service  facility  has  uni ini  ted 
waiting  space. 

Assumption  3.  There  can  be  any  number  of  servers  at  each 
node.  They  are  identical  independent  servers,  each  serving 
one  customer  at  a  tiae. 

Assumption  4.  Customers  are  selected  for  service  at  each 
facility  according  to  the  first  come,  first-served 
discipline. 

Assumption  5.  There  can  be  any  number  of  customer 
classes ,  but  customers  cannot  change  classes.  Moreover, 
much  of  the  analysis  in  QNA  is  done  for  the  aggregate  or 
typical  customer. 

Assumption  6.  Customers  can  be  created  or  combined  at 
the  nodes,  e.g.  an  arrival  can  cause  more  than  one 
departure. 


(Whitt,  1983a:2781-2782) 


QNA  uses  two  parameters  to  characterize  the  arrival  process  and  the 
service  times — one  to  describe  the  rate  and  the  other  to  describe  the 
variability  (Whitt,  1983a:2782). 

Required  Inputs.  QNA  allows  several  different  formats  for 
entering  the  necessary  information.  In  general,  the  Information  that 
must  be  supplied  is  as  follows:  (1)  the  number  of  nodes  in  the  network, 
(2)  the  number  of  servers  at  each  node,  (3)  the  external  arrival  rate  to 
each  node,  (4)  the  variability  parameter  of  the  external  arrival  process 
to  each  node,  (5)  the  mean  service  time  at  each  node,  (6)  the  squared 
coefficient  of  variation  of  the  service-time  distribution  at  each  node, 
and  (7)  the  Markovian  routing  of  entities  within  the  network. 

QNA  Outputs.  QNA  will  provide  the  steady-state  congestion 
measures  for  each  node  in  the  network.  The  main  congestion  measure  is 
the  mean  waiting  time  (before  beginning  service),  but  QNA  also  generates 
an  entire  probability  distribution  for  the  waiting  time.  QNA  will  also 
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provide  the  probability  that  the  server  is  busy  at  an  arbitrary  tine, 
the  expected  nuaber  of  entities  in  the  node,  the  probability  that  an 
entity  is  delayed,  and  the  conditional  delay  given  that  the  server  is 
busy  (Whitt,  1983a:2802-2807) . 

QNA  will  also  calculate  the  approximate  congestion  aeasures  for 
the  network  as  a  whole.  It  provides  congestion  aeasures  representing 
the  systea  view  (e.g.  throughput  and  nuaber  of  entities  in  the  network) 
and  congestion  measures  representing  the  customer  view  (e.g.  nuaber  of 
nodes  visited  and  response  tines)  (Whitt,  1983a:2807). 

Performance  of  QUA.  Whitt  (1983b)  describes  the  performance 
of  QNA  and  compares  the  congestion  measures  to  those  obtained  through 
simulation  and  the  standard  Markovian  algorithm  (which  is  represented  by 
the  M/M/m  equations  given  in  Chapter  I).  He  tested  the  performance  of 
QNA  on  a  variety  of  queueing  networks  from  a  single  GI/C/1  queue  to  a 
packet-switched  communication-network  model.  The  results  of  Whitt's 
study  demonstrated  the  importance  of  the  variability  parameter  used  in 
QNA  to  estimate  the  congestion  measures  of  networks  that  do  not  satisfy 
the  assumptions  of  the  Jackson  network.  Furthermore,  when  the  Jackson 
network  assumptions  are  satisfied,  then  the  approximations  used  in  QNA 
yield  the  exact  measures. 
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III .  Methodology 

The  objective  of  this  study  was  to  compare  several  control 
variates  that  have  shown  proaising  results  for  queueing  network 
siaulation.  This  was  done  with  the  ultiaate  goal  in  aind  of  providing 
the  siaulation  coaaunity  with  soae  guidance  for  selecting  control 
variates  that  will  lead  to  significant  reductions  in  the  variance  of  the 
estiaated  responses  while  not  introducing  bias  to  the  estimates. 

The  major  obstacle  to  achieving  the  stated  goal  was  to  design  an 
experiment  that  is  general  enough  so  that  the  results  are  applicable  to 
a  wide  range  of  queueing  network  models.  And,  at  the  sane  tine,  the 
size  of  the  experiment  aust  be  manageable  so  that  it  can  be  completed  in 
the  alloted  time. 

Description  of  the  Queueing  Network 

The  first  step  in  keeping  with  the  above  considerations  was  to 
select  a  queueing  network.  The  author  decided  to  select  a  single 
network  that  is  small  in  terms  of  the  number  of  nodes  but  complex  enough 
to  incorporate  many  aspects  of  queueing  networks  in  general. 

The  basic  form  of  the  network  was  based  upon  the  third  network 
that  was  used  by  Sharon  (1986).  It  is  an  open  network  which  consists  of 
four  nodes,  each  of  which  has  a  single  server.  Entities  arrive  from 
outside  the  network  to  the  first  service  center  (or  node)  according  to  a 
Poisson  process  with  arrival  rate  A  =  1.  At  each  service  center,  the 
queue  capacity  is  infinite,  and  the  queue  discipline  is  FIFO.  From  the 
first  node  entities  can  go  to  either  node  2  or  node  3  with  given 
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probabilities.  Entities  from  nodes  2  and  3  proceed  to  node  4.  Finally, 
entities  leaving  node  4  may  loop  back  to  the  first  node  or  may  exit  the 
system.  The  basic  structure  of  the  queueing  network  is  illustrated 
below  in  Figure  3.1. 


Furtheraore,  the  service  tiaes  at  the  first  three  nodes  were 
distributed  exponentially.  This  decision  was  made  to  siaplify  the 
experlaent,  rather  than  introduce  another  variable  to  the  experimental 
design.  In  effect,  this  does  not  limit  the  applicability  of  the  results 
because  the  superposition  of  independent  arrival  processes  approaches  a 
Poisson  process  in  the  liait  no  matter  how  they  were  originally 
distributed. 

However,  the  service  tine  distribution  at  the  fourth  node  was 
varied  to  obtain  the  desired  variability  in  the  service  tiaes,  since  an 
exponential  distribution  is  Halted  to  a  coefficient  of  variation  equal 
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to  one.  The  coefficient  of  variation  is  a  Measure  of  the  variability  of 
a  distribution  and  is  given  by  the  square  root  of  the  variance  divided 
by  the  Bean.  The  naxinua  entropy  distribution  was  chosen  to  generate 
service  tines  with  a  coefficient  of  variation  less  than  one  and  the 
hyperexponential  distribution  was  chosen  to  generate  service  tines  with 
a  coefficient  of  variation  greater  than  one.  There  is  aore  detail  about 
these  distributions  later  in  this  chapter. 

The  Response  Variables 

Next,  two  response  variables  were  selected  for  the  analysis.  One 
of  then  was  the  average  sojourn  tine  for  an  entity  to  pass  through  the 
network,  which  can  be  approxinated  by  finding  the  solution  to  the 
approxiaating  Jackson  Network.  It  can  also  be  solved  by  finding  the 
approxinate  solution  to  the  exact  network  using  QNA. 

The  other  response  variable  selected  was  a  quantile  at  the  fourth 
node  representing  the  probability  that  the  nuaber  in  queue  exceeds  soae 
threshhold  value.  The  threshhold  value  chosen  was  twice  the  aean  nuaber 
in  queue  at  equilibriun.  For  practical  purposes,  QNA  was  used  to  find 
the  threshhold  value  by  taking  the  next  highest  integer  of  the  following 
result:  2  X  (EN  -  p )  ,  where  EN  is  the  expected  nuaber  at  the  fourth 
node,  and  p  is  the  associated  traffic  intensity. 

Quantiles  are  not  easily  estiaated,  even  though  knowledge  of  then 
nay  prove  to  be  iaportant.  To  illustrate,  consider  a  coanunication 
network  where  the  queues  represent  buffers  for  incoaing  Messages.  And, 
suppose  that  one  proposes  to  deteraine  the  capacity  of  the  buffers  in 
this  network  by  finding  the  aean  nuaber  in  queue  at  equilibriun  for  a 
systea  with  infinite  capacity  and  siaply  doubling  the  result.  The  user 
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of  the  network  >ay  then  want  to  know  how  often  the  capacity  of  the 
buffers  are  exceeded  to  deteraine  whether  or  not  the  buffers'  capacities 
should  be  increased. 


The  Control  Variables 

Next,  examples  of  both  types  of  control  variates,  internal  and 
external,  were  selected.  The  author  picked  two  different  types  of 
internal  control  variates  (standardized  routing  controls  and 
standardized  work  variables)  and  two  different  types  of  external  control 
variates  (analytic  Jackson  controls  and  analytic  approximations) . 

The  Internal  Control  Variates.  Two  standardized  routing  controls, 
one  for  each  of  the  two  probabilistic  branchings  in  the  network,  were 
selected  for  this  research.  For  further  reference,  let  R,x  denote  the 
routing  control  for  the  proportion  of  entities  that  took  the  path  from 
node  1  to  node  3  (as  opposed  to  node  1  to  node  2),  and  let  R«,  denote 
the  routing  control  for  the  proportion  of  entities  that  took  the  path 
from  node  4  back  to  node  1.  These  standardized  routing  controls  are 
given  by  Eq  (2.37). 

In  addition  to  the  two  routing  controls,  four  standardized  work 
variables,  one  for  each  of  the  service  nodes,  were  selected.  The  four 
work  variables,  denoted  by  W,,  W»,  W,,  and  W,,  are  given  by  Eq  (2.36). 

The  External  Control  Variates.  Also,  two  external  control 
variates  under  each  of  the  two  different  approaches  alluded  to  earlier 
were  selcted  for  this  research.  These  external  controls  were  the 
steady-state  expected  sojourn  time  for  the  entire  network  and  the 
steady-state  expected  waiting  time  in  the  fourth  queue. 
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The  analytic  Jackson  control  variates  and  the  analytic 
approximations  were  found  using  Bell  Laboratories'  Queueing  Network 
Analyzer  (QNA).  The  following  input  parameters  to  the  simulation  model 
were  used  to  generate  the  "known"  means  of  the  external  control 
variates:  the  external  arrival  rate,  the  mean  service  times  at  each 
node,  the  squared  coefficient  of  variation  of  the  service  times  at  each 
node,  and  the  probabilistic  routing  matrix.  Under  the  assumption  that 
the  squared  coefficient  of  variation  of  the  service  times  are  all  equal 
to  one  yields  the  M/M/1  (or  Jackson  Network)  results.  On  the  other 
hand,  using  the  input  squared  coefficient  of  variation  for  the  service- 
times  at  the  fourth  node  (the  only  one  that  violates  the  above 
asumptlon)  yields  the  C/C/1  (or  approximate)  results. 

The  observed  values  of  the  external  control  variates  were  found 
using  the  observed  values  of  these  same  parameters  as  inputs  to  QNA. 

The  input  files  to  QNA  for  generating  the  "known"  means  at  the  first 
design  point  are  provided  in  Appendix  B,  and  the  "known"  means  of  the 
external  control  variates  at  all  design  points  are  summarized  in  a  table 
in  Appendix  C. 

Selecting  a  Service-Time  Distribution 

In  almost  all  of  the  literature  examined,  there  has  not  been  any 
rationale  specified  for  selecting  a  particular  distributional  form  for 
the  interarrival  times  or  service  times.  Rather,  most  researchers 
select  a  few  of  the  more  commonly  used  distributions  and  include  the 
type  of  distribution  as  a  variable  in  their  research.  However,  in  this 
case,  in  keeping  with  the  goal  of  maintaining  enough  generality,  the 
selection  of  a  specific  distribution  might  bias  the  results.  Therefore, 


the  naxiaua  entropy  distribution  was  selected  for  this  research  because 
it  is  the  least  biased. 

However,  because  of  the  difficulties  associated  with  generating 
randoa  nuabers  according  to  the  aaxlaua  entropy  distribution,  the  author 
decided  to  use  it  for  generating  service  tines  at  the  fourth  service 
center  only.  As  indicated  earlier,  the  fourth  service  center  is  of 
particular  interest  because  of  the  selection  of  the  quantile 
representing  the  probability  that  the  nuaber  in  queue  at  the  fourth 
service  center  exceeds  twice  the  expected  nuaber  under  steady-state 
conditions. 

Unfortunately,  the  author  discovered  that  the  naxiaua  entropy 
distribution  could  not  be  used  to  generate  a  nonnegative  random  variate 
with  a  coefficient  of  variation  (c.)  greater  than  one.  (The  author  did 
not  explore  whether  or  not  the  aaxiaua  entropy  distribution  can  have  a 
coefficient  of  variation  greater  than  one  over  the  real  nunbers.) 
Instead,  the  hyperexponential  distribution  was  selected  for  such  cases. 
This  selection  is  also  suggested  from  graphical  representations  of  the 
distributions.  At  the  lower  setting  of  the  coefficient  of  variation 
(0.5)  the  aaxiaua  entropy  distribution  looks  like  a  normal  distribution 
that  was  truncated  at  the  origin  (negative  tines  are  impossible).  As 
the  coefficient  of  variation  increases  to  one  the  graph  looks 
exponential.  In  fact,  when  the  coefficient  of  variation  equals  one  the 
naxiaua  entropy  distribution  reduces  to  the  exponential  distribution. 

The  hyperexponential  distribution  with  coefficient  of  variation  greater 
than  one  takes  on  an  exponential  shape  but  with  a  longer  tail. 


Therefore,  the  aaxiaua  entropy  distribution  was  used  to  generate 
randoa  variates  with  a  coefficient  of  variation  less  than  one  and  the 
hyperexponential  distribution  was  used  to  generate  randoa  variates  with 
a  coefficient  of  variation  greater  than  one.  The  functional  foras  of 
the  density  functions  for  these  two  distributions  are  given  below: 

Maxiaua  Entropy:  f(x)  =  exp(-l  -  A,  -  A,x  -  A*x*)  (3.1) 

Hyperexponential:  f(x)  »  c,/(3,  *exp(-x/e, )  +  Ct/^-expt-x/e,)  (3.2) 

Generating  Randoa  Variates  Osina  Selected  Distributions 

The  author  used  the  siaulation  prograaaing  language  SLAM  II  to 
code  up  the  network  aodel.  SLAM  II  has  built-in  functions  to  generate 
randoa  variates  according  to  several  distributions.  However,  it  does 
not  have  any  built-in  functions  for  generating  the  aaxiaun  entropy  or 
hyper exponential  distributions.  Therefore,  the  author  had  to  build 
FORTRAN  subroutines  using  well-known  techniques  to  generate  the  randoa 
variates  desired.  The  SLAM  II  Network  code  and  FORTRAN  subroutines  are 
provided  for  the  reader  in  Appendix  A. 

Maxiaua  Entropy  Distribution.  The  deveiopaent  of  the  parameters 
(A,,  A,,  and  A,)  for  the  aaxiaua  entropy  distribution  given  the  first 
two  aoaents  (|i,  and  |i2)  is  no  trivial  aatter.  All  of  the  autoaated 
search  procedures  tried  have  not  been  successful,  but  a  prograa  using 
MINOS  is  currently  being  pursued.  However,  the  author  was  able  to  find 
paraaeters  through  a  soaewhat  manual  search  process  using  a  routine  to 
nuaerically  evaluate  the  appropriate  integrals  to  reproduce  the  desired 
aoaents  to  within  4  parts  in  10,000.  Given  the  density  function  f(x)  as 
defined  in  Eq  (3.1)  the  three  integrals  are  defined  below: 
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(3.3) 


r 


f(x)  dx  =  1 


r 


x-f(x)  dx  =  J1, 


x)  dx  =  »1, 


(definition  of  a  density  function) 


(definition  of  the  first  moment) 


(definition  of  the  second  moment) 


(3.4) 


(3.5) 


However,  having  found  the  parameters  still  leaves  the  problem  of 
generating  random  variates  according  to  the  maximum  entropy  distribution 
with  those  parameters.  To  generate  the  random  variates  according  to  the 
desired  maximum  entropy  distribution  the  author  decided  to  use  the 
acceptance-rejection  method  (Law  and  Kelton,  1982:250-252). 

The  efficiency  of  the  acceptance-rejection  method  is  determined  by 
the  area  between  the  density  fuctions  of  the  majorizing  distribution  and 
the  original  distribution  of  interest.  As  the  area  decreases  the 
efficiency  of  the  method  increases.  If  f(x)  and  g(x)  are  density 
functions  of  two  different  distributions,  g(x)  majorizes  f(x)  if  and 
only  if  g(x)  >  f(x)  at  every  x  where  f(x)  is  defined. 

Using  HathCAD  the  author  was  able  to  visually  fit  a  majorizing 
distribution  to  each  maximum  entropy  distribution.  The  majorizing 
distributions  were  developed  using  the  density  function  of  a  Normal 
distribution  that  was  truncated  at  the  origin.  To  insure  that  they 
truly  majorized  the  maximum  entropy  distribution  they  were  checked  by 
creating  a  table  of  values  in  NathCAD  and  adjusting  them  as  necessary. 
Figures  3.2  through  3.5  are  graphical  representations  of  the  density 
functions  for  the  four  maximum  entropy  distributions  that  were  used  in 
this  research  with  their  respective  majorizing  distributions. 
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Figure  3.2.  Service-Time  Distribution  at  Node  #4 
with  m,  =  0.45  and  (i,  =  0.253125 


Figure  3.3.  Service-Time  Distribution  at  Node  #4 
with  ii,  =  0.81  and  ji,  =  0.820125 
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Figure  3.4.  Service-Tine  Distribution  at  Node  #4 
with  n,  a  0.375  and  p,  ■  0.1757812 


Figure  3.5.  Service-Tine  Distribution  at  Mode  #4 
with  |i,  *  0.675  and  |it  -  0.5695312 


Hvoer exponential  Distribution.  Given  the  first  two  moments,  the 


parameters  of  the  hyperexponent lal  distribution  can  be  solved  for 
directly  fron  the  four  equations  given  below. 


c,8,  +  cf8*  *  U, 

(first  moment) 

(3.6) 

2c, 8, 1  *  2c,8, 1 

«  |i,  (second  moment) 

(3.7) 

c,  +  c,  =  1 

(3.8) 

c,8,  *  c*8i 

(3.9) 

Eqs  (3.6)  and  (3.7)  were  derived  fron  the  definitions  of  the  first  and 
second  moments,  respectively.  Eq  (3.8)  was  derived  fron  the  fact  that 
the  integral  of  any  density  function  equals  one.  And,  Eq  (3.9)  is  added 
to  provide  a  unique  solution.  Simultaneously  solving  Eqs  (3.6)  through 
(3.9)  in  terns  of  the  first  two  nonents  leads  to  the  following  results: 


c. 

*  (0.25  -  O.Sm.Vj*,)''*  +  0.5 

(3.10) 

c. 

*  1  -  c, 

(3.11) 

e, 

■  0.5u,/c, 

(3.12) 

e. 

■  0.5|l,/Ca 

(3.13) 

A  hyper exponential  distribution  that  satisifies  Eq  (3.9)  is  said  to  be 
balanced.  The  reader  is  referred  to  pages  139-147  of  Kleinrock  (1975) 
for  more  information. 

To  generate  random  variates  according  to  the  hyperexponential 
distribution  the  author  selected  the  composition  method  (Law  and  Kelton, 
1982:247-249).  Figures  3.6  through  3.9  are  graphical  representations  of 
the  density  functions  of  the  four  hyperexponential  distributions  that 
were  used  in  this  research. 
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Lag-1  Correlation  o£  the  Interdeparture  Times 

Prior  to  any  experimentation  with  the  network  described  earlier, 
the  author  experimented  with  a  single-server  queue  to  test  for  the 
strength  of  the  lag-1  correlation  of  the  interdeparture  times  of  a  C/C/1 
queue.  The  queue  capacity  was  infinite  and  the  queue  discipline  was 
FIFO.  Furthermore,  the  maximum  entropy  and  hyperexponential 
distributions  were  used  to  generate  the  Interarrival  times  and  the 
service  times  of  the  entities.  The  author  assumed  a  first-order 
autoregressive  time  series  model  for  the  autocorrelation  function  of  the 
departure  process. 

Experimental  Design.  Also,  the  author  used  a  2*  factorial  design 
to  explore  the  effect  of  the  traffic  intensity  and  the  variability  of 
the  interarrival-tine  and  service-tine  distributions  on  the  first-order 
autoregressive  parameter.  The  levels  of  the  three  factors  (traffic 
intensity,  coefficient  of  variation  of  the  interarrival  tines,  and 
coefficient  of  variation  of  the  service  times)  in  the  experiment  are 
summarized  in  the  table  below. 

Table  3.1.  Levels  of  Factors  in  2*  Experiment 
of  a  G/G/l  Queue 


Low  Value 

High  Value 

Parameter 

(-1) 

(♦1) 

P 

0.5 

0.9 

c. 

O.S 

2.5 

c. 

0.5 

2.5 

3.U 


To  eliminate  the  Initialization  bias,  the  author  discarded  the 
first  five  thousand  observations  of  the  interdeparture  tines.  This 
decision  was  nade  based  upon  a  visual  interpretation  of  a  tine- 
persistent  plot  of  the  average  waiting  tine  in  the  queue.  Then,  the 
next  one  thousand  observations  of  the  Interdeparture  tines  were  used  to 
fit  a  first-order  autoregressive  tine-series  nodel. 

Exnerinental  Results.  Using  PROC  ARINA  under  SAS  the  resulting 
first-order  autoregressive  paraneter  ranged  fron  -0.117  to  +0.078. 
Although  these  extrene  values  proved  to  be  statistically  significant 
(nost  of  the  intervening  values  did  not),  they  are  not  practically 
significant. 

This  result  is  significant  to  the  developnent  of  analytic 
approxinations  for  estinating  the  perfornance  neasures  of  a  network. 

That  is,  since  the  arrival  process  to  a  queue  in  a  network  is  the 
superposition  of  an  external  arrival  process  and  the  departure  processes 
of  other  queues  in  the  network,  then  this  results  validates  the 
assumption  that  the  interarrival  tlaes  are  independent.  Otherwise,  the 
approxlaation  formulas  would  have  to  incorporate  the  dependency  between 
the  interarrival  tines  of  the  entities  to  the  queues  in  the  network. 
Furthermore,  the  author's  findings  provide  the  reason  why  Whitt  (1984) 
did  not  find  any  laproveaent  when  he  incorporated  the  lag-1  correlation 
of  the  interdeparture  tines  in  his  two-aoaent  approxinations. 

Experimental  Design 

The  experinental  design  for  conparing  the  control  variates  over 
the  network  Illustrated  in  Figure  3.1  was  a  2*  factorial  design.  The 
four  factors  selected  were  the  traffic  intensity  of  the  network,  the 
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coefficient  of  variation  of  the  service  tinea  at  the  fourth  node,  and 
the  two  routing  variables  (r,,  and  r«,).  The  levels  of  these  four 
factors  are  given  in  the  table  below. 

Table  3.2.  Levels  of  Factors  in  2*Experiaent 
of  the  Queueing  Network 


Parameter 

Low  Value 
(-1) 

High  Value 
(+1) 

P 

0.5 

0.9 

c. 

0.5 

2.5 

r.i 

0.2 

0.4 

0.1 

0.25 

Given  the  external  arrival  rate  to  the  first  node  (which  was  set 
equal  to  one)  and  the  routing  matrix,  the  effective  arrival  rates  to 
each  of  the  four  nodes  in  the  network  can  be  solved  for  from  the  traffic 
rate  equations  given  by 


n 

A,  .  A,,  +  l  A,r, ,  (3.14) 

i»l 


where 

A,  *  the  effective  arrival  rate  to  node  j 

A#J  *  the  external  arrival  rate  to  node  j 

A,  «  the  effective  arrival  rate  to  node  i 

r , ,  *  the  routing  probability  from  node  1  to  node  j 

n  >  the  nuaber  of  nodes  in  the  network 


Next  the  aean  service  tine  (r)  for  each  node  (j)  can  be  solved  for  by 


Finally,  Che  first  two  moments  of  the  service-time  distribution  at  the 
fourth  node  can  be  solved  for  by 


l»,  *  t*  (3.16) 

P.  -  p,*(l  +  c.*)  (3.17) 


The  values  of  the  input  parameters  to  the  simulation  model  at  each  of 
the  design  points  are  summarized  in  Table  3.3. 


Table  3.3.  Values  of  the  Input  Parameters 
to  the  Queueing  Network  Simulation  Model 


Design  — Factor  Settings—  - Mean  Service  Times 


Point 

P 

c. 

r,  * 

r4i 

f. 

ft 

1* 

*.(p. ) 

Pi 

1 

0.5 

0.5 

0.4 

0.10 

0.450 

0.75000 

1.1250 

0.450 

0.2531250 

2 

0.9 

0.5 

0.4 

0.10 

0.810 

1.35000 

2.0250 

0.810 

0.8201250 

3 

0.5 

2.5 

0.4 

0.10 

0.450 

0.75000 

1.1250 

0.450 

1.4681250 

4 

0.9 

2.5 

0.4 

0.10 

0.810 

1.35000 

2.0250 

0.810 

4.7567250 

5 

0.5 

0.5 

0.2 

0.10 

0.450 

0.56250 

2.2500 

0.450 

0.2531250 

6 

0.9 

0.5 

0.2 

0.10 

0.810 

1.01250 

4.0500 

0.810 

0.8201250 

7 

0.5 

2.5 

0.2 

0.10 

0.450 

0.56250 

2.2500 

0.450 

1.4681250 

8 

0.9 

2.5 

0.2 

0.10 

0.810 

1.01250 

4.0500 

0.810 

4.7567250 

9 

0.5 

0.5 

0.4 

0.25 

0.375 

0.67500 

0.9375 

0.375 

0.1757812 

10 

0.9 

0.5 

0.4 

0.25 

0.675 

1.12500 

1.6875 

0.675 

0.5695312 

11 

0.5 

2.5 

0.4 

0.25 

0.375 

0.67500 

0.9375 

0.375 

1.0195313 

12 

0.9 

2.5 

0.4 

0.25 

0.675 

1.12500 

1.6875 

0.675 

3.3032813 

13 

0.5 

0.S 

0.2 

0.25 

0.375 

0.46875 

1.8750 

0.375 

0  1757812 

14 

0.9 

0.5 

0.2 

0.25 

0.675 

0.84375 

3.3750 

0.675 

0.5695312 

15 

0.5 

2.5 

0.2 

0.25 

0.375 

0.46875 

1.8750 

0.375 

1.0195313 

16 

0.9 

2.5 

0.2 

0.25 

0.675 

0.84375 

3.3750 

0.675 

3.3032813 

Note  that  at  the  design  points  where  c.  *  0.5-  the  service-time 
distribution  at  the  fourth  node  was  the  maximum  entropy  distribution; 
and,  at  those  design  points  where  c,  *  2.5,  the  service-time 
distribution  at  the  fourth  node  was  the  hyperexponential  distribution. 


The  parameters  for  generating  the  hyperexponential  distribution  were 
derived  in  a  FORTRAN  subroutine  in  the  model  using  Eqs  (3.10)  through 
(3.13).  But,  the  parameters  for  generating  the  maximum  entropy 
distribution  were  derived  outside  of  the  model.  These  parameters  were 
then  input  to  the  model  and  are  given  In  Table  3.4  below. 

Table  3.4.  Parameters  for  Generating  Random  Variates 
According  to  the  Maximum  Entropy  Distribution 


Parameters  of  the  Density  Functions 


Maximum  Entropy  Dist 

Majorizing 

Dist 

»** 

^8 

x, 

c 

<r* 

0.2531250 

-0.028146 

-7.000 

8.1980 

1.12 

m 

Finn 

0.8201250 

0.557130 

-3.882 

2.5265 

1.06 

EKH1 

IK 

0.375 

0.1757812 

-0.210450 

-8.400 

11.8050 

1.13 

0 . 350 

o.os 

0.5695312 

0.375020 

-4.659 

3.6385 

1.10 

0.640 

0.15 

Collecting  Data 

Now  that  we  have  selected  a  queueing  network,  the  response 
variables,  the  control  variables,  the  interarrival-time  and  service-time 
distributions,  and  an  experimental  design,  we  need  to  decide  how  much 
data  to  collect  and  how  to  generate  the  data.  The  first  question,  “how 
much  data  to  collect,”  is  related  to  estimating  the  bias  introduced  to 
the  controlled  estimators  by  forming  confidence  intervals  about  them  and 
counting  what  percentage  of  the  confidence  intervals  cover  the  "known” 
value  of  the  response.  The  second  question,  "how  to  generate  the  data,” 
refers  to  choosing  a  technique  that  will  yield  unbiased  independent 
observations  of  the  uncontrolled  estimators. 
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Determining  How  Much  Data  to  Collect.  The  author  decided  to  use 
twenty  observations  (that  is  twenty  runs  of  the  simulation  model)  to 
obtain  each  estimate  of  the  mean  controlled  responses.  Then,  about  the 
estimated  controlled  mean  response  a  single  95X  confidence  interval  was 
formed.  Now,  to  obtain  a  reasonable  estimate  of  the  actual  coverage 
percentages  of  the  confidence  interval  the  author  decided  that  this 
required  a  minimum  of  two  hundred  confidence  intervals.  Therefore,  four 
thousand  runs  (twenty  runs  per  confidence  interval  X  two  hundred 
confidence  intervals)  were  required  to  determine  the  bias  introduced  to 
the  controlled  responses  at  each  of  the  sixteen  design  points. 

Determining  How  tg.  Generate  the  Data.  Knowing  that  we  need  four 
thousand  runs  of  the  simulation  model  at  each  of  the  sixteen  design 
points  poses  some  problems  in  terms  of  the  computer  time  required  to 
make  all  these  runs.  The  two  major  considerations  here  are  eliminating 
the  initialization  bias  and  deciding  whether  to  use  the  method  of 
independent  replications  or  batch  means;  and,  both  of  these 
considerations  are  related. 

Independent  Replications  vs  Batch  Means .  The  method  of 
Independent  replications  will  yield  Independent  observations  since  each 
run  of  the  simulation  model  uses  different  random  numbers.  However,  the 
transient  period  of  the  simulation  must  be  discarded  for  each  run.  On 
the  other  hand,  the  transient  period  is  discarded  only  once  for  the 
batch  means  approach,  but  the  observations  may  be  autocorrelated.  The 
practitioner  must  then  determine  the  appropriate  batch  size  necessary  to 
pass  some  statistical  tests  of  independence.  Any  good  text  on 
simulation  modeling  should  have  more  information  about  the  batch  means 
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approach;  a  few  of  then  are  Kleijnen  (1974),  Law  and  Kelton  (1982),  and 
Pri taker  (1986). 

In  general,  when  one  is  faced  with  Baking  so  many  runs  of  a 
simulation  model  the  batch  neans  approach  is  preferred  to  independent 
replications  because  of  the  savings  in  computer  time.  But,  the  effect 
of  the  batch  means  approach  upon  the  correlation  between  the  response 
variable(s)  and  the  control  variable(s)  has  not  been  studied  (or,  at 
least,  the  author  has  not  found  any  literature  reporting  such  a  study). 
In  many  situations  there  is  a  lag  time  in  the  correlation  between  the 
response  variable(s)  and  the  control  variable(s). 

For  example,  suppose  for  our  queueing  network  we  select  the  mean 
sojourn  time  for  the  response  variable  and  the  mean  service  time  at  the 
first  node  for  the  control  variable.  Then,  when  the  observed  average 
service  time  at  the  first  node  is  high  compared  to  its  known  mean,  we 
would  expect  the  average  sojourn  time  to  be  higher  than  its  true  mean. 
But,  we  would  also  expect  that  the  average  sojourn  time  is  not  affected 
iaediately  by  some  longer  than  average  service  times  at  the  first  node. 
Therefore,  when  using  the  batch  means  approach,  if  the  batch  size  is  not 
large  enough  to  capture  the  correlation  between  the  response  variable! s) 
and  the  control  variable(s)  we  may  not  achieve  a  variance  reduction,  or 
bias  may  be  introduced  to  the  controlled  estimates. 

For  these  reasons  the  author  decided  to  use  independent 
replications  instead  of  the  batch  means  approach. 

Eliminating  the  Initialization  Bias.  Since  we  are 
interested  in  steady-state  results,  the  transient  (or  warm-up)  period  of 
the  simulation  must  be  discarded.  The  author  made  some  trial  runs  of 
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the  simulation  and  plotted  a  time-persistent  average  of  the  sojourn 
time.  The  author  made  a  visual  determination  that  the  transient  period 
was  essentially  over  after  about  three  thousand  time  units  when  the 
simulation  was  started  at  empty  and  idle  conditions.  But,  this  would  be 
too  costly  in  terms  of  computer  time  to  discard  this  much  data  for  each 
run  of  the  model.  Therefore,  the  author  decided  to  start  the  simulation 
model  essentially  at  steady  state  by  initializing  the  number  of  entities 
in  each  queue  to  the  expected  number  as  given  by  QNA.  Then,  in  order  to 
provide  some  randomization,  the  observations  generated  during  the  first 
one  hundred  time  units  were  discarded. 

Weighting  the  Observations 

Since  observations  were  collected  on  so  many  different  variables 
(two  response  variables,  six  internal  control  variables,  and  four 
external  control  variables)  it  was  more  convenient  to  make  each 
simulation  run  of  the  same  fixed  length  (which  was  one  thousand  time 
units)  rather  than  try  to  obtain  the  same  number  of  observations  for 
each  variable.  That  is,  each  observation  simulation  run  of  the 
variables  under  study  (resulting  from  one  run  of  the  simulation)  are 
averages  of  the  individual  observations  within  a  simulation  run  over  one 
thousand  time  units;  and,  although  the  time  interval  is  constant,  the 
number  of  observations  during  the  fixed  time  interval  (within  a 
simulation)  varies  from  one  run  to  the  next.  Therefore,  to  obtain 
unbiased  estimates  of  discrete  performance  variables  (such  as  the  mean 
sojourn  time)  we  need  to  weight  the  observations  from  each  simulation 
run  accordingly. 
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For  example,  each  run  of  the  simulation  yields  one  value  of  the 
average  sojourn  tiae,  which  itself  is  an  average  of  all  the  individual 
sojourn  times  of  the  entities  completing  the  network  during  the  fixed 
time  interval  of  one  thousand  time  units.  Then,  the  uncontrolled 
estimate  of  the  mean  sojourn  time  is  a  weighted  average  of  twenty 
averages,  and  the  weights  were  derived  from  the  number  of  individual 
►  observations  that  comprised  each  of  the  twenty  averages  respectively. 

l 

The  weighted  average  of  any  random  variable  X  is  given  by 

K  K 

X  -  £  «.X./  2  w«  (3.18) 

I  •  I  ft  •  I 

where  K  is  the  number  of  observations  (simulation  runs)  and  w,  is  the 
weighting  coefficient  on  the  ith  observation  of  X.  The  weighting 
coefficients  are  given  by 


*»« 
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where  n,  is  the  number  of  individual  observations  (within  a  simulation) 
that  make  up  the  ith  observation  of  X.  Similarly,  the  weighted  sample 
variance  of  X  is  given  by 


K 

I 

*1 

s*  = 

l  l  w»X, *  - 
1  •  1 

2  «»x< 

I  •  s 

/[X(K-1)1 

(3.20) 

Statistics  Used  to.  Compare  Control  Variates 

Finally,  we  need  to  specify  the  measures  to  be  used  to  compare  the 


results  of  the  different  control  variables  on  the  response  variables. 


First  of  all,  the  control  variates  were  compared  on  the  basis  of 
efficiency  (or  the  size  of  the  variance  reduction)  as  given  by  the 
variance  ratio  (VR).  Recall  fro*  Chapter  II  that  the  VR  is  the  ratio  of 
the  variance  of  the  controlled  estimator  when  the  optimal  control 
coefficient  is  unknown  (and  must  be  estimated)  to  the  variance  of  the 
uncontrolled  estimator. 

Secondly,  the  potential  bias  introduced  to  the  controlled 
estimates  was  measured  by  the  percentage  of  the  total  number  of 
confidence  intervals  about  the  controlled  estimates  that  covered  the 
grand  mean  of  the  4000  uncontrolled  observations. 

The  Experimental  Procedure 

In  this  section,  the  author  traces  his  steps  in  running  the 
programs  and  generating  the  data.  Recall  that  four  thousand  runs  of  the 
simulation  model  were  made  at  each  of  the  sixteen  design  points  and 
that  each  run  of  the  simulation  model  was  eleven  hundred  time  units  in 
length. 

Because  of  CPU  time  limits,  one  thousand  runs  of  the  simulation 
model  were  submitted  at  a  time;  therefore,  four  submissions  were 
required  to  complete  a  single  design  point.  Each  of  these  job 
submissions  used  appro,  imately  one  hour  of  CPU  time  on  a  VAX  86S0 
computer  and  generated  five  output  files.  The  four  replications  of 
these  output  files  for  each  design  point  were  then  appended  to  one 
another. 

The  contents  of  these  five  output  files  are  as  follows:  (1) 
JACKSON. IN  contained  the  records  of  the  parameters  necessary  for  input 
to  QNA  to  solve  for  the  observed  values  of  the  two  external  control 
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variates  under  the  assuaption  that  all  the  service  variability 
paraaeters  equalled  one;  (2)  QNA. IN  was  siailar  to  JACKSON. IN  except 
that  the  observed  values  of  the  service  variability  paraaeters  were 
included;  (3)  RESPONSE. OUT  contained  the  nuaber  of  entities  coapleting 
the  network  and  the  two  response  variables;  (4)  ROUTING. OUT  contained 
the  two  standardized  routing  control  variables  with  their  respective 
nuaber  of  observations;  and  (5)  WORK. OUT  contained  the  total  nuaber  of 
service  coapletions  and  the  four  standardized  work  variables.  The 
reader  aay  exaaine  the  FORTRAN  code  that  generated  these  files  in 
Appendix  A. 

Next,  the  author  used  a  slightly  aodified  version  of  the  QNA 
software  was  to  produce  the  two  output  files  JACKSON.OUT  and  QNA. OUT 
which  contained  the  observations  on  the  external  control  variates 
generated  froa  their  respective  input  files.  (The  author's 
aodif ications  to  QNA  were  to  suppress  the  noraal  detailed  output  and  to 
report  the  values  of  the  two  variables  of  interest  only). 

The  five  files  with  the  ".OUT"  extension  (along  with  the  two  files 
M16.0UT  and  C16.0UT  which  contained  the  "known”  aeans  of  the  external 
control  variates)  then  becaae  the  input  files  to  the  prograa  called 
CONTROL.  The  files  M16.0UT  and  C16.0UT  are  given  in  table  foraat  in 
Appendix  C,  and  the  FORTRAN  code  for  the  prograa  CONTROL  is  provided  in 
Appendix  A.  The  prograa  CONTROL  was  written  to  calculate  the 
uncontrolled  <.3  well  as  the  controlled  estiaates  of  the  aeans  of  the  two 
response  variables  against  each  of  the  ten  control  variates.  It  also 
provided  the  variance  ratios  and  95X  confidence  liaits  about  the 
controlled  responses.  Each  output  record  was  based  on  twenty  records  of 
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input  data  (representing  twenty  runs  of  the  simulation)  and  from  here  is 
referred  to  as  a  nacro-replication;  and,  there  were  two  hundred  macro- 
replications  at  each  design  point.  The  output  data  (macro-replications) 
were  stored  in  four  files:  (1)  SOJOURN. VR  which  contained  the 
uncontrolled  estimates,  the  controlled  estimates,  and  the  variance 
ratios  for  the  response  variable  sojourn  time;  (2)  SOJOURN. Cl  which 
contained  the  95X  confidence  limits  for  the  response  variable  sojourn 
time;  (3)  QUANTILE. VR  which  contained  the  uncontrolled  estimates,  the 
controlled  estimates,  and  the  variance  ratios  for  the  response  variable 
representing  the  probability  that  the  number  in  queue  exceeded  twice  the 
expected  number;  and  (4)  QUANTILE. Cl  which  contained  the  respective  95X 
confidence  limits. 

Finally,  the  two  SOJOURN  files  and  the  two  QUANTILE  files  were 
separately  input  to  the  program  RESULTS  which  calculated  the  minimum, 
mean,  and  maximum  values  of  the  uncontrolled  estimates,  the  controlled 
estimates,  and  the  variance  ratios  over  the  two  hundred  macro- 
replications.  RESULTS  also  counted  the  number  of  times  the  grand 
uncontrolled  estimate  of  the  mean  response  fell  within  the  95X 
confidence  limits  of  the  two  hundred  macro-replications  and  reports  the 
percentage  that  do  so  as  an  estimate  of  the  actual  coverage  of  the 
confidence  intervals.  The  output  files  produced  by  the  program  RESULTS 
are  provided  in  Appendix  D.  Also,  the  FORTRAN  code  for  the  program 
RESULTS  is  provided  in  Appendix  A. 
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IV.  Results 


The  results  of  the  experimentation  with  control  variates  on  the 
open  queueing  network  described  in  Chapter  III  are  given  in  this 
chapter.  The  variance  ratios  and  confidence  interval  coverages  are 
reported  for  ten  control  variates  against  two  response  variables.  The 
output  files  from  the  program  RESULTS  are  provided  in  Appendix  D. 

Table  4.1  lists  the  grand  means  of  the  uncontrolled  simulation 
response  of  sojourn  time  and  the  analytic  solutions  provided  by  the 
M/M/1  and  C/C/1  formulas.  Note  that  the  analytic  solutions  were  used  as 
the  "known"  means  of  the  external  control  variate  of  sojourn  time  as 
given  in  Appendix  C. 

Since  the  M/M/1  analytic  results  assume  that  the  squared 
coefficient  of  variation  of  the  service  times  at  the  fourth  node  is 
equal  to  one,  we  would  expect  a  difference  between  the  results  obtained 
with  the  M/M/1  formulas  and  the  results  obtained  from  the  simulation 
model.  However,  it  is  interesting  to  note  the  significant  differences 
between  the  simulation  results  and  those  obtained  using  the  C/C/1 
approximations.  These  differences  are  given  in  terms  of  percentages  in 
Table  4.1  below. 

The  differences  between  the  estimated  mean  sojourn  times  for  an 
entity  to  complete  the  network  were  more  pronounced  as  the  traffic 
intensity  of  the  network  increased;  and  likewise,  these  differences 
increased  as  the  coefficient  of  variation  of  the  service  times  at  the 
fourth  node  increased.  These  findings  are  summarized  in  Table  4.2 
below. 
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Table  4.1.  Analytic  Results  vs.  Simulation  Results 
For  Estimating  Mean  Sojourn  Time  in  Hetwork 


- Analytic 

M/M/1 

Results - 

C/C/1 

Simulation 

Results 

X  Delta 

1 

0.400000E+01 

0.38U77Ef01 

0. 381331E+01 

0.04 

2 

0 . 360000E+02 

0. 329330E+02 

0.292945E+02 

-11.05 

3 

0.400000E+01 

0.531762E+01 

0.521644E+01 

-1.90 

4 

0. 360000E+02 

0.574688E+02 

0.457255E+02 

-20.43 

5 

0.400000E+01 

0. 381173E+01 

0 . 301449E+O1 

0.07 

6 

0. 360000E+02 

0. 329329E+02 

0.287763E+02 

-12.62 

7 

0.400000E+01 

0.531792E+01 

0 .521470E+01 

-1.94 

8 

0 . 360000E+02 

0.574698E+02 

0.451113E+02 

-21.50 

9 

0.417391E+01 

0.398269E+01 

0 . 398874E+01 

0.15 

10 

0 . 360000E+02 

0. 327806E+02 

0.299427E+02 

-8.66 

11 

0.417391E+01 

0.551237E+01 

0 .530441E+01 

-3.77 

12 

0 . 360000E+02 

0 .  585327E+02 

0.470301E+02 

-19.65 

13 

0.400000E+01 

0.380878E+01 

0.381095E+01 

0.06 

14 

0 . 360000E+C2 

0 . 327797E+02 

0 . 292389E+02 

-10.80 

IS 

0.400000E+01 

0.S33846E+01 

0.512254E+01 

-4.04 

16 

0. 360000E+02 

0.585391E+02 

0 . 456687E+02 

-21.99 

*X  Delta  is  the  percent  difference  between  the  C/C/1  results  and 
the  simulation  results. 


Table  4.2.  Effect  of  p  and  c,  on  X  Delta 
Between  Analytic  C/C/1  and  Simulation  Results 


Average 

P 

c. 

X  Delta 

0.5 

0.5 

SKVVB 

2.5 

-2.91 

0.9 

ImEM 

-10.78 

0.9 

2.5 

-20.89 

Because  the  differences  increase  as  the  traffic  intensity  of  the 
network  increases  and  as  the  variability  of  the  service  times  at  the 
fourth  node  increases,  this  may  indicate  one  of  three  things,  or 
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possibly  a  combination  of  the  three.  First,  the  C/C/1  approximations 
used  in  QNA  may  be  less  accurate  as  the  traffic  intensity  increases  and 
as  the  coefficient  of  variation  of  the  service  times  increases.  Or, 
there  may  be  significant  initialization  bias  in  the  simulation  results 
at  the  higher  traffic  intensity  and  the  higher  variability  of  the 
service  times.  Finally,  the  fixed  time  interval  of  the  simulation  may 
have  been  too  short  to  provide  good  estimates  at  the  higher  traffic 
intensity  and  higher  variability  of  the  service  times. 

The  author’s  first  suspicion  that  the  C/C/1  approximations  used  in 
QNA  might  be  a  significant  contributor  to  these  differences  seems  to  be 
justified  in  light  of  the  findings  of  Uhitt  (1983b).  As  to  the 
initialization  bias,  the  author  started  the  simulation  with  the  expected 
number  in  each  queue.  However,  these  numbers  were  provided  by  QNA;  and 
therefore,  under  the  assumption  that  QNA's  results  were  less  accurate  as 
the  traffic  intensity  and  service  variability  increased  an  initial  bias 
was  probably  Introduced.  Although  further  experimentation  with  the 
simulation  model  is  required  to  prove  any  of  these  suspicions,  knowing 
these  differences  exist  proved  helpful  in  interpreting  the  resulting 
variance  ratios  and  coverages  of  the  confidence  intervals. 

The  average  variance  ratios  achieved  using  each  control  variate 
against  the  first  response  variable,  sojourn  time,  are  given  in  Table 
4.3.  At  the  bottom  of  each  column  in  Table  4.3  appears  the  grand 
average  variance  ratio  for  each  control  variate  across  the  sixteen 
design  points.  Recall  that  the  smaller  the  variance  ratio,  the  greater 
the  variance  reduction.  Table  4.4  provides  the  average  variance  ratios 
when  the  first  two  factors  (traffic  intensity  and  service  variability) 
are  held  constant. 
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Table  4.3.  Variance  Ratios  Achieved 
With  Sojourn  Tine  as  Response  Variable 


- Internal  Controls -  - External  Controls - 

—Routing . . Work  Variables -  — M/M/1 -  — G/G/l - 


Pt 

R 1 1 

R<  1 

w, 

wt 

w. 

w. 

SOJ 

w. 

SOJ 

w. 

1 

0.937 

0.868 

0.955 

0.949 

0.949 

0.955 

0.418 

0.768 

0.385 

0.751 

2 

0.942 

0.907 

0.944 

0.946 

0.9S1 

0.944 

0.697 

0.645 

0.698 

0.642 

3 

0.941 

0.925 

0.956 

0.951 

0.954 

0.956 

0.468 

0.386 

0.251 

0.317 

4 

0.942 

0.933 

0.947 

0.942 

0.943 

0.948 

0.662 

0.615 

0.601 

0.611 

5 

0.870 

0.886 

0.957 

0.957 

0.954 

0.957 

0.378 

0.797 

0.355 

0.783 

6 

0.941 

0.917 

0.960 

0.947 

0.951 

0.960 

0.771 

0.650 

0.772 

0.644 

7 

0.928 

0.914 

0.950 

0.951 

0.949 

0.950 

0.508 

0.422 

0.279 

0.369 

8 

0.940 

0.946 

0.951 

0.948 

0.939 

0.952 

0.704 

0.620 

0.609 

0.617 

9 

0.950 

0.762 

0.939 

0.948 

0.943 

0.939 

0.355 

0.662 

0.334 

0.643 

10 

0.951 

0.859 

0.943 

0.951 

0.948 

0.943 

0.708 

0.604 

0.710 

0.602 

11 

0.946 

0.874 

0.947 

0.949 

0.945 

0.947 

0.450 

0.416 

0.266 

0.387 

12 

0.950 

0.917 

0.949 

0.953 

0.946 

0.949 

0.689 

0.672 

0.671 

0.673 

13 

0.878 

0.777 

0.946 

0.949 

0.942 

0.946 

0.331 

0.712 

0.315 

0.700 

14 

0.929 

0.858 

0.952 

0.945 

0.948 

0.952 

0.748 

0.633 

0.750 

0.627 

15 

0.935 

0.899 

0.939 

0.958 

0.948 

0.939 

0.448 

0.404 

0.253 

0.379 

16 

0.942 

0.929 

0.945 

0.944 

0.951 

0.945 

0.726 

0.624 

0.655 

0.624 

Avg  0.933 

0.886 

0.949 

0.949 

0.948 

0.949 

0.566 

0.602 

0.494 

0.586 

Table  4.4.  Average  Variance  Ratios  for  the  First  Response 
Across  Traffic  Intensity  and  Service  Variability 


Control 

-Variance 

Ratios  at 

Factor  Levels 

(/*,  c.)  — 

Variate 

(0.5,0. 5) 

(0.9, 0.5) 

(0.5, 2. 5) 

(0.9, 2. 5) 

Ri  * 

0.909 

0.941 

0.938 

0.943 

R«  i 

0.823 

0.885 

0.903 

0.931 

W, 

0.949 

0.950 

0.948 

0.948 

W, 

0.951 

0.947 

0.952 

0.947 

W. 

0.947 

0.949 

0.949 

0.945 

w« 

0.949 

0.950 

0.948 

0.948 

S(M/M/1 ) 

0.371 

0.731 

0.469 

0.695 

W1M/M/1) 

0.735 

0.633 

0.407 

0.633 

S(G/G/1) 

0.347 

0.733 

0.262 

0.634 

W( G/G/l) 

0.719 

0.629 

0.363 

0.631 

The  results  listed  in  Table  4.3  show  that  among  the  internal 
control  variates  R«t  achieved  the  smallest  variance  ratio  and  the  others 
were  roughly  equal.  However,  the  external  control  variates  achieved 
significantly  greater  variance  ratios  that  any  individual  Internal 
control  variate.  Also,  among  the  external  control  variates,  those  based 
upon  the  G/G/l  formulas  achieved  smaller  variance  ratios  than  those 
based  upon  the  M/M/1  formulas. 

Furthermore,  one  can  see  from  Table  4.4  that  except  for  R,,  the 
average  variance  ratios  of  the  internal  control  variates  are  not  greatly 
affected  by  changes  in  the  traffic  intensity  or  service  variability. 
However,  the  average  variance  ratios  for  the  external  control  variates 
are  greatly  affected  by  both  the  traffic  intensity  and  service 
variability. 

In  a  similar  manner,  the  average  variance  ratios  achieved  using 
each  control  variate  against  the  second  response  variable,  the  quantile 
representing  the  probability  of  the  number  in  the  fourth  queue  exceeding 
twice  the  expected  number,  are  given  in  Table  4.5.  Likewise,  Table  4.6 
provides  the  average  variance  ratios  when  the  first  two  factors  (traffic 
intensity  and  service  variability)  are  held  constant. 

The  same  general  observations  hold  true  for  the  results  of  the 
variance  ratios  achieved  against  the  second  response  variable  as  for  the 
first  response  variable.  Additionally,  in  comparing  the  variance  ratios 
between  the  two  response  variables,  the  variance  ratios  are  smaller  for 
the  first  response  variable. 

Next,  Tables  4.7  through  4.10  are  the  respective  analogs  to  Tables 
4.3  through  4.6  with  the  confidence  interval  coverages  as  the  statistic 
for  comparison  in  place  of  the  variance  ratios. 


Table  4.5.  Variance  Ratios  Achieved 
With  Quantile  at  4th  Node  as  Response  Variable 


- Internal  Controls -  - External  Controls - 

— Routing —  - Work  Variables -  — M/M/1 -  — G/G/l - 


Pt 

R|  3 

t 

w. 

w. 

w. 

w« 

SOJ 

w. 

SOJ 

w. 

1 

0.939 

0.909 

0.950 

0.946 

0.952 

0.951 

0.705 

0.443 

0.702 

0.421 

2 

0.945 

0.918 

0.946 

0.944 

0.948 

0.946 

0.918 

0.723 

0.920 

0.721 

3 

0.946 

0.943 

0.952 

0.951 

0.951 

0.952 

0.580 

0.268 

0.263 

0.233 

4 

0.951 

0.956 

0.945 

0.947 

0.941 

0.945 

0.837 

0.752 

0. 775 

0.748 

5 

0.947 

0.904 

0.950 

0.953 

0.951 

0.950 

0.779 

0.446 

0.797 

0.418 

6 

0.951 

0.939 

0.953 

0.950 

0.948 

0.953 

0.940 

0.73S 

0.941 

0.732 

7 

0.943 

0.935 

0.945 

0.948 

0.952 

0.945 

0.669 

0.270 

0.310 

0.243 

8 

0.944 

0.934 

0.957 

0.951 

0.948 

0.957 

0.864 

0.747 

0.775 

0.746 

9 

0.943 

0.838 

0.948 

0.952 

0.942 

0.948 

0.612 

0.396 

0.624 

0.374 

10 

0.953 

0.914 

0.951 

0.9S5 

0.953 

0.951 

0.900 

0.713 

0.901 

0.711 

11 

0.947 

0.934 

0.945 

0.946 

0.946 

0.945 

0.611 

0.271 

0.272 

0.239 

12 

0.946 

0.935 

0.948 

0.950 

0.943 

0.948 

0.849 

0.743 

0.797 

0.740 

13 

0.948 

0.840 

0.946 

0.944 

0.947 

0.946 

0.691 

0.382 

0.712 

0.359 

14 

0.946 

0.905 

0.9S3 

0.949 

0.956 

0.953 

0.923 

0.704 

0.924 

0.700 

15 

0.947 

0.948 

0.950 

0.955 

0.947 

0.950 

0.653 

0.256 

0.300 

0.238 

16 

0.947 

0.949 

0.952 

0.947 

0.948 

0.952 

0.852 

0.712 

0.773 

0.707 

Avg  0.946 

0.919 

0.949 

0.949 

0.948 

0.949 

0.774 

0.535 

0.674 

0.521 

Table  4.6.  Average  Variance  Ratios  for  the  Second  Response 
Across  Traffic  Intensity  and  Service  Variability 


Control 

-Variance 

Ratios  at 

Factor  Levels 

(P,  c.)~ 

Variate 

(0.5, 0.5) 

(0.9, 0.5) 

(0.5, 2. 5) 

(0.9, 2. 5) 

ft.. 

0.944 

0.949 

0.946 

0.947 

1 

0.873 

0.919 

0.940 

0.944 

H, 

0.949 

0.951 

0.948 

0.951 

w. 

0.949 

0.949 

0.950 

0.949 

w. 

0.948 

0.951 

0.949 

0.945 

w. 

0.949 

0.951 

0.948 

0.951 

S (M/M/1) 

0.697 

0.920 

0.628 

0.850 

W(M/M/1 ) 

0.417 

0.719 

0.266 

0.739 

S (G/G/l) 

0.709 

0.922 

0.286 

0.780 

W(C/G/1) 

0.393 

0.716 

0.238 

0.735 
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Table  4.7.  Coverage  Percentages  of  95X  Confidence  Interval 
About  Controlled  Estimate  of  Mean  Sojourn  Time 


- Internal 

Controls - 

- External 

Controls - 

—Routing - 

— 

-Work  Variables- 

— M/M/1 - 

— C/G/l - 

Pt 

ft.. 

R#  i 

w, 

w. 

w. 

w. 

SOJ 

w. 

SOJ 

w. 

1 

0.930 

0.885 

0.920 

0.925 

0.935 

0.920 

0.550 

0.830 

0.545 

0.810 

2 

0.935 

0.955 

0.940 

0.935 

0.930 

0.940 

0.800 

0.695 

0.800 

0.700 

3 

0.945 

0.945 

0.950 

0.950 

0.945 

0.950 

0.570 

0.510 

0.365 

0.430 

4 

0.930 

0.935 

0.910 

0.900 

0.915 

0.910 

0.615 

0.640 

0.580 

0.625 

5 

0.910 

0.925 

0.935 

0.960 

0.955 

0.935 

0.495 

0.815 

0.475 

0.820 

6 

0.920 

0.925 

0.945 

0.940 

0.935 

0.945 

0.810 

0.745 

0.810 

0.720 

7 

0.920 

0.895 

0.930 

0.925 

0.905 

0.925 

0.610 

0.545 

0.420 

0.510 

8 

0.945 

0.955 

0.945 

0.930 

0.950 

0.945 

0.655 

0.690 

0.590 

0.695 

9 

0.915 

0.765 

0.895 

0.930 

0.905 

0.895 

0.465 

0.695 

0.440 

0.710 

10 

0.910 

0.895 

0.890 

0.925 

0.895 

0.890 

0.785 

0.685 

0.780 

0.660 

11 

0.935 

0.905 

0.930 

0.935 

0.935 

0.930 

0.580 

0.580 

0.465 

0.525 

12 

0.915 

0.905 

0.925 

0.935 

0.920 

0.925 

0.680 

0 . 725 

0.685 

0.720 

13 

0.880 

0.800 

0.895 

0.910 

0.915 

0.895 

0.430 

0.765 

0.435 

0.740 

14 

0.945 

0.900 

0.945 

0.930 

0.935 

0.945 

0.815 

0.745 

0.815 

0.745 

15 

0.930 

0.925 

0.915 

0.925 

0.925 

0.915 

0.565 

0.515 

0.385 

0.495 

16 

0.915 

0.905 

0.895 

0.915 

0.905 

0.890 

0.655 

0.650 

0.630 

0.665 

Avg  0.924 

0.901 

0.923 

0.929 

0.925 

0.922 

0.630 

0.677 

0.576 

0.661 

Table  4.8.  Average  Coverages  for  the  First  Response 
Across  Traffic  Intensity  and  Service  Variabiliy 


Control 

Variate 

-Variance 
(0.5, 0.5) 

Ratios  at  Factor  Levels  (p,c.)  — 
(0.9, 0.5)  (0.5, 2. 5)  (0.9, 2. 5) 

R., 

0.909 

0.928 

0.933 

0.926 

R«  i 

0.844 

0.919 

0.917 

0.925 

Ws 

0.911 

0.930 

0.931 

0.919 

w. 

0.931 

0.933 

0.934 

0.920 

W, 

0.928 

0.924 

0.927 

0.923 

W, 

0.911 

0.930 

0.930 

0.918 

S(M/M/1 ) 

0.485 

0.803 

0.581 

0.651 

W1M/N/1) 

0.776 

0.717 

0.538 

0.676 

S( C/G/l) 

0.474 

0.801 

0.409 

0.621 

W(G/G/1 ) 

0.770 

0.706 

0.490 

0.676 
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Table  4.9.  Coverage  Percentages  of  95X  Confidence  Interval 
About  Controlled  Estimate  of  Fourth  Mode  Quantile 


- Internal  Controls -  - External  Controls - 

— Routing —  - Work  Variables -  — M/M/1 -  — G/G/l - 


Pt 

K  1  3 

ft#  a 

w, 

w. 

Wj 

w« 

SOJ 

w. 

SOJ 

w. 

1 

0.905 

0.890 

0.900 

0.910 

0.895 

0.900 

0.750 

0.570 

0.770 

0.560 

2 

0.940 

0.930 

0.930 

0.920 

0.935 

0.930 

0.915 

0.765 

0.915 

0.760 

3 

0.940 

0.930 

0.955 

0.940 

0.930 

0.955 

0.660 

0.435 

0.365 

0.375 

4 

0.850 

0.865 

0.855 

0.840 

0.845 

0.855 

0.735 

0.685 

0.68S 

0.675 

5 

0.910 

0.885 

0.910 

0.905 

0.905 

0.910 

0.780 

0.545 

0.795 

0.510 

6 

0.925 

0.935 

0.930 

0.925 

0.930 

0.930 

0.915 

0.800 

0.915 

0.800 

7 

0.940 

0.930 

0.925 

0.930 

0.915 

0.925 

0.750 

0.485 

0.485 

0.435 

8 

0.870 

0.865 

0.865 

0.860 

0.865 

0.865 

0.755 

0.720 

0.670 

0.715 

9 

0.930 

0.850 

0.910 

0.940 

0.925 

0.910 

0.690 

0.495 

0.690 

0.465 

10 

0.935 

0.920 

0.930 

0.940 

0.925 

0.930 

0.930 

0.760 

0.930 

0.750 

11 

0.965 

0.960 

0.935 

0.955 

0.930 

0.935 

0.710 

0.470 

0.485 

0.430 

12 

0.87S 

0.860 

0.880 

0.885 

0.875 

0.880 

0.755 

0.730 

0.745 

0.730 

13 

0.915 

0.830 

0.905 

0.900 

0.905 

0.905 

0.705 

0.510 

0.715 

0.505 

14 

0.960 

0.925 

0.935 

0.920 

0.940 

0.935 

0.915 

0.765 

0.915 

0.760 

15 

0.935 

0.950 

0.905 

0.925 

0.945 

0.905 

0.730 

0.470 

0.500 

0.445 

16 

0.850 

0.850 

0.855 

0.850 

0.830 

0.855 

0.705 

0.670 

0.655 

0.655 

Avg 

0.915 

0.898 

0.908 

0.909 

0.906 

0.908 

0.775 

0.617 

0.702 

0.598 

Table  4.10.  Average  Coverages  for  the  Second  Response 
Across  Traffic  Intensity  and  Service  Variabiliy 


Control 

Variate 

-Variance 
(0.5,0. 5) 

Ratios  at  Factor  Levels  (p, c.)-- 
(0.9, 0.5)  (0.5, 2. 5)  (0.9, 2. 5) 

R,, 

0.915 

0.940 

0.945 

0.861 

R«, 

0.864 

0.928 

0.943 

0.860 

w, 

0.906 

0.931 

0.930 

0.864 

W, 

0.914 

0.926 

0.938 

0.859 

W, 

0.907 

0.933 

0.930 

0.854 

w. 

0.906 

0.931 

0.930 

0.864 

S(H/M/1) 

0.731 

0.919 

0.713 

0.737 

W(M/M/1 ) 

0.530 

0.773 

0.465 

0.701 

S(G/G/1 ) 

0.743 

0.919 

0.459 

0.689 

W(G/G/1 ) 

0.510 

0.767 

0.421 

0.694 
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A  general  stateaent  that  can  be  Bade  about  the  coverage 
percentages  of  the  95%  confidence  intervals  about  the  controlled  mean 
responses  is  that  the  coverage  worsens  as  the  variance  ratio  decreases. 
This  result  is  not  satisfactory,  because  we  would  like  to  be  able  to 
achieve  large  variance  reductions  (i.e.  sea 11  variance  ratios)  and  not 
bias  the  controlled  responses  (i.e.  good  coverages  of  the  confidence 
intervals) . 

Finally,  another  way  of  exaeining  the  saee  results  is  presented  in 
Tables  4.11  through  4.14.  Here,  these  four  tables  give  the  average 
effects  for  the  four  factors.  An  effect,  froa  regression  analysis,  is 
siaply  the  average  value  of  the  response  when  a  factor  is  at  its  high 
setting  ainus  the  average  value  of  the  response  when  the  saae  factor  is 
at  its  low  setting.  The  way  to  interpret  these  tables  is  that  the 
factors  with  the  larger  effects  in  aagnitude  have  a  greater  iapact  on 
the  statistic  (variance  ratios  or  coverages  of  confidence  intervals). 

Froa  Tables  4.11  through  4.14  one  can  see  that  the  four  factors 
have  a  relatively  saall  effect  upon  the  results  using  the  internal 
control  variates;  however,  two  factors  (the  traffic  intensity  and  the 
service  variability)  have  a  relatively  large  effect  upon  the  results 
using  the  external  control  variates. 
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Table  4.11.  Average  Effects  of  the  Four  Factors 
on  the  Variance  Ratios  of  the  First  Response 


- Effects  of  the  Four  Factors - 

Control  Traffic  C(  of  4th  r,,  r« 


Variate 

Intensity 

Server 

Prob. 

Prob. 

^  t  $ 

0.009 

0.008 

0.012 

0.003 

1 

0.023 

0.031 

-0.005 

IPBK  ^  i 

w, 

0.000 

-0.001 

-0.001 

w. 

-0.002 

0.000 

-0.001 

W, 

0.000 

-0.001 

0.000 

»4 

0.000 

-0.001 

-0.001 

WEL  W  i : 

S  <  M/M/l ) 

0.147 

0.016 

-0.010 

W( M/M/1) 

0.031 

-0.082 

-0.006 

S<C/C/1) 

0.189 

-0.046 

-0.005 

W(G/C/1) 

0.044 

-0.088 

-0.007 

Table  4.12.  Average  Effects  of  the  Four  Factors 
on  the  Variance  Ratios  of  the  Second  Response 


- Effects  of  the  Four  Factors - 

Control  Traffic  Ct  of  4th  r,,  r« 


Variate 

Intensity 

Server 

Prob. 

Prob. 

0.001 

0.000 

0.000 

0.001 

R«. 

0.012 

0.023 

0.000 

-0.011 

w, 

0.001 

0.000 

-0.001 

0.000 

w. 

0.000 

0.000 

0.000 

0.001 

w. 

0.000 

-0.001 

-0.001 

-0.001 

w. 

0.001 

0.000 

-0.001 

0.000 

S< M/M/l) 

0.111 

-0.035 

-0.022 

-0.013 

W(M/M/1 ) 

0.194 

-0.033 

0.004 

-0.013 

S(G/G/1 ) 

0.177 

-0.141 

-0.017 

-0.011 

W(G/C/1 ) 

0.205 

-0.034 

0.003 

-0.012 
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Table  4.13.  Average  Effects  of  the  Four  Factors 
on  the  Coverages  of  the  First  Response 

- Effects  of  the  Four  Factors - 

Control  Traffic  C(  of  4th  r,,  r«, 


Variate 

Intensity 

Server 

Prob. 

Prob. 

R  i  j 

0.003 

0.006 

0.003 

-0.006 

R«, 

0.021 

0.020 

-0.003 

-0.026 

w, 

0.002 

0.002 

-0.003 

-0.012 

W, 

-0.003 

-0.003 

0.000 

-0.004 

W, 

-0.002 

0.000 

-0.003 

-0.008 

w. 

0.002 

0.002 

-0.002 

-0.012 

S( M/M/1) 

0.097 

-0.014 

0.001 

-0.008 

W(M/M/1 ) 

0.020 

-0.070 

-0.007 

-0.007 

S(C/C/1) 

0.135 

-0.061 

0.006 

0.003 

W(G/G/1 ) 

0.031 

-0.077 

-0.013 

-0.003 

Table  4.14.  Average  Effects  of  the  Four  Factors 
on  the  Coverages  of  the  Second  Response 


Effects  of  the  Four  Factors 


Control 

Traffic 

C,  of  4th 

r,. 

r., 

Variate 

Intensity 

Server 

Prob. 

Prob. 

R,, 

-0.015 

-0.012 

0.002 

0.005 

R., 

-0.005 

0.003 

0.002 

-0.005 

w, 

-0.010 

-0.011 

0.004 

-0.001 

w. 

-0.017 

-0.011 

0.007 

0.005 

w. 

-0.013 

-0.014 

0.002 

0.003 

w. 

-0.010 

-0.011 

0.004 

-0.001 

S(M/M/1 ) 

0.053 

-0.050 

-0.007 

-0.008 

W(M/M/1 ) 

0.120 

-0.034 

-0.003 

-0.008 

S(G/G/1 ) 

0.102 

-0.128 

-0.004 

0.002 

W(G/G/1 ) 

0.132 

-0.041 

-0.005 

-0.006 
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V.  Conclusions 


The  purpose  of  this  study  was  to  compare  several  control  variates 
for  queueing  network  simulation.  The  author's  goal  was  to  provide  the 
simulation  community  with  some  guidance  for  selecting  control  variates 
that  will  lead  to  significant  reductions  in  the  variance  of  the 
estimated  responses  that  do  not  introduce  significant  bias.  The 
conclusions  of  this  research  are  important  because  they  add  to  the  body 
of  knowledge  a  simulation  practitioner  can  draw  upon  when  seiecting  a 
variance  reduction  technique.  Also,  these  results  indicate  that  further 
research  in  this  area  is  warranted. 

The  results  of  this  research  are  consistent  with  previous  findings 
demonstrating  the  potential  variance  reduction  that  can  be  obtained  when 
using  control  variates.  However,  this  research  went  beyond  most  of  the 
previous  efforts  in  measuring  the  bias  introduced  to  the  controlled 
estimates  by  way  of  the  percentage  of  confidence  intervals  about  the 
controlled  responses  that  covered  the  grand  averages  of  the  uncontrolled 
responses. 

Another  novel  addition  to  this  research  was  the  use  of  the  maximum 
entropy  distribution  to  generate  service  times,  which  is  the  least- 
biased  distribution  when  only  the  first  two  moments  of  the  distribution 
are  known.  However,  finding  the  parameters  to  the  maximum  entropy 
distribution  was  more  difficult  than  expected,  and  generating  random 
variates  according  to  the  maximum  entropy  distribution  was  an  involved 
process.  The  author  recommends  that  future  research  efforts  look  at 
comparing  the  Eriang  distribution  (which  is  much  easier  to  generate)  to 
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the  maximum  entropy  distribution  for  generating  random  variates  with  a 
coefficient  of  variation  less  than  one.  The  hyperexponential 
distribution,  however,  is  the  recommended  distribution  for  generating 
random  variates  with  a  coefficient  of  variation  greater  than  one  when 
only  the  first  two  moments  of  the  distribution  are  known. 

The  use  of  the  Queueing  Network  Analyzer  (QNA)  to  provide  external 
control  variates  was  also  a  new  approach.  However,  the  results  showed 
that  QNA  could  yield  estimates  of  the  congestion  measures  as  much  as 
twenty  percent  in  error  to  the  simulation  estimates,  which  might  have 
been  a  potential  source  of  bias.  Future  research  is  also  needed  to 
determine  exactly  how  these  differences  can  be  reduced. 

The  way  to  reduce  the  errors  in  the  estimates  of  the  congestion 
measures  of  QNA  is  to  use  better  approximations.  One  suggestion  is  to 
use  the  approximation  for  the  expected  waiting  time  in  the  GI/G/m  queue 
given  by  Kimura  (1986),  since  he  reported  achieving  better  results  as 
compared  to  those  achieved  using  QNA  (the  reader  i3  referred  to  Chapter 
II).  The  author  pursued  this  suggestion  by  computing  the  expected 
waiting  times  at  the  fourth  node  over  the  sixteen  design  points  using 
Kimura's  formula.  However,  the  largest  percentage  change  from  QNA's 
result  was  only  -1.29X;  therefore,  using  Kimura's  formula  in  place  of 
the  one  in  QNA  would  not  have  helped  in  this  case.  These  results  are 
presented  for  the  reader  in  Table  5.1. 

In  general,  the  external  control  variates  achieved  smaller 
variance  ratios  than  the  internal  control  variates;  however,  the 
coverages  of  the  confidence  intervals  about  the  controlled  responses 
were  worse.  The  range  of  the  average  variance  ratios  for  the  internal 
control  variates  was  0.886  to  0.949  with  coverages  from  0.898  to  0.929 
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Table  S.l.  Waiting  Tines  at  Fourth  Node 
Computed  Using  QNA's  and  Kimura's  Approximations 


DESIGN 

—EXPECTED  WAITING  TIMES— 

%  CHANGE 

POINT 

QNA'S 

KIMURA'S 

FROM  QNA 

1 

0.281197E+00 

0.281135E+00 

0.022 

2 

0.45558SE+01 

0.455531E+01 

0.012 

3 

0 . 163162E+01 

0.163416E+01 

-0.1S6 

4 

0 . 264290E+02 

0 . 264305E+02 

-0.006 

5 

0.281160E+00 

0.281053E+00 

0.038 

6 

0 . 455572E+01 

0 . 455566E+01 

0.001 

7 

0 . 163188E+01 

0 . 163623E+01 

-0.267 

8 

0 . 264300E+02 

0. 264320E+02 

-0.008 

9 

0.234173E+00 

0.233935E+00 

0.102 

10 

0 . 379484E+01 

0 . 379464E+01 

0.005 

11 

0. 136079E+01 

0. 137057E+01 

-0.719 

12 

0.220362E+02 

0 . 220441E+02 

-0.036 

13 

0 . 234015E+00 

0 . 233590E+00 

0.182 

14 

0 . 379416E+01 

0 . 379390E+01 

0.007 

15 

0 . 136189E+01 

0. 137946E+01 

-1.290 

16 

0 . 220409Ef02 

0. 220515E+-02 

-0.048 

for  the  95%  confidence  intervals.  The  range  of  the  average  variance 
ratios  for  the  external  control  variates  was  0.494  to  0.774  with 
coverages  from  0.576  to  0.775  for  the  95%  confidence  intervals. 

More  specifically,  smaller  variance  ratios  were  achieved  against 
the  first  response  variable,  the  average  sojourn  time  for  the  network, 
than  for  the  second  response  variable,  the  probability  of  the  number  in 
the  fourth  queue  exceeding  twice  the  expected  number.  This  was  expected 
since  one  would  expect  the  control  variates  chosen  to  have  a  greater 
correlation  with  the  first  response  variable  in  general. 

Also,  the  standardized  routing  controls  achieved  smaller  variance 
ratios  than  the  standardized  work  variables.  Particularly,  any  route 
which  feeds  entities  back  through  a  portion  of  the  network  (thereby 
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increasing  the  congestion)  is  a  good  candidate  for  an  internal  control 
variable. 

When  comparing  the  results  of  the  internal  control  variates  to  the 
external  control  variates  one  should  keep  in  Bind  that  the  standardized 
internal  controls  are  generally  Independent  of  one  another.  Therefore, 
aultiple  internal  controls  could  be  used  to  achieve  even  larger  variance 
reductions  provided  the  nunber  of  replications  is  sufficient  to  overcome 
the  loss  factor.  It  is  also  inportant  to  recognize  that  the  external 
control  variates  tend  to  introduce  bias  to  the  controlled  estimates. 
Another  area  for  future  research  would  be  into  reducing  this  bias  since 
the  potential  for  variance  reductions  is  so  promising.  Some  potential 
ways  to  reduce  the  bias  are  to  use  jacknife  estimators  or  some  other 
estimators  that  do  not  shrink  the  confidence  interval  width  as  much. 

Furthermore,  future  efforts  should  look  into  the  effect  of  the 
batch  means  approach  on  the  control  variate  technique.  The  batch  means 
approach  provides  a  means  of  reducing  the  cost  of  obtaining  several 
thousand  observations  when  trying  to  estimate  the  coverage  percentages 
of  the  confidence  intervals.  However,  exactly  how  the  batch  size 
affects  the  correlation  between  the  response  variable(s)  and  the  control 
variable(s)  is  not  known;  and  therefore,  the  author  used  the  independent 
replications  approach. 

An  important  observation  to  note  is  that  the  standardized  internal 
controls  were  less  influenced  by  the  changes  in  the  factor  settings  than 
the  external  controls,  which  in  turn  yielded  more  consistent  results 
over  the  experimental  design.  In  particular,  the  internal  controls  are 
robust  for  the  traffic  intensity. 


The  author  recommends  using  internal  control  variates  since  they 
demonstrated  better  coverage;  however,  multiple  controls  will  be 
necessary  to  achieve  a  significant  reduction.  The  use  of  external 
controls  is  very  promising  since  a  single  control  can  achieve  a  very 
significant  variance  reduction;  however,  further  research  should  be  done 
to  determine  ways  to  reduce  the  bias  they  Introduce  to  the  estimated 
responses. 

When  selecting  internal  control  variates  for  use  in  a  queueing 
network  simulation  the  author  recommends  the  standardized  routing 
controls  over  the  standardized  work  variables,  especially  when  there  are 
probabilistic  branches  in  the  network  that  have  a  significant  effect  on 
the  response(s).  For  example,  in  the  author's  research  the  entities 
completing  the  fourth  node  could  either  exit  the  system  or  be  fed  back 
to  the  first  node.  Any  deviation  from  the  expected  number  of  entities 
being  fed  back  has  a  significant  impact  upon  the  estimate  of  the  average 
sojourn  time  in  the  network;  and  therefore,  the  routing  control  R,, 
achieved  a  significantly  greater  variance  reduction  than  the  other 
internal  control  variates. 

On  a  final  note,  when  comparing  the  control  variates  the  smaller 
variance  ratios  indicate  a  larger  variance  reduction.  The  increased 
number  of  runs  that  would  be  required  to  achieve  the  same  variance 
reduction  without  the  control  variate(s)  can  be  calculated  by  rounding 
the  result  given  by  Eq  (5.1)  to  the  next  highest  integer. 

K...  iTi.m  =  K-tCl/VR)  -  1]  (5.1) 
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For  example,  If  we  made  twenty  runs  of  the  simulation  and  achieved  a 
variance  ratio  of  0.90,  then  we  would  need  3  additional  runs.  However, 
If  we  achieved  a  variance  ratio  of  0.50,  then  we  would  need  20 
additional  runs  (for  a  total  of  40  runs)  to  achieve  the  same  amount  of 
variance  reduction  without  control  variates.  Thus,  significant  savings 
in  computer  time  can  be  realized  through  the  use  of  control  variates. 
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Appendix  A:  Computer  Source  Code 


SLAW  II  Network  Code 


GEN, QUEUEING  NETWORK, JOHN  TOMICK, 10/20/88, 1000, N,N,Y/Y,N,N, 72; 
LIMITS, 4, 2, 500; 

INTLC,XX(1)*0.5,XX(2)*0.5,XX(3)=0.45,XX(4)=0.253125; 
INTLC,XX(6)*0 .45 ,XX( 7 )*0. 75 ,XX(8)»i . 1 25 , XX  <  9 ) = 1 00 ; 

INTLC,XX( 10)*-0. 028146, XX( 11 )*-7.0,XX(12)=*8. 198, XX(13)*1. 12; 
INTLC,XX<14)*0.425,XX(15>*0.07,XX(20)*Q.4,XX(2i)*0.1,XX(22)=l. ; 


NETWORK; 

CREATE, EXPON<l.,l), 0,1; 

EVENT, 1; 

QUE1  ASS IGN , ATRI B<  2 ) = EX PON! XX( 6 ) , 2 ) ; 
EVENT, 2; 

QUEUE! 1 ) , 1 ; 

ACTIVITY/1 ,ATRIB( 2) ; 

EVENT, 3; 

GOON  1  * 

ACTIVITY,, XX(20),QUE3; 
ACTIVITY; 

QUE2  ASSIGN, ATRI B( 2 )=>EXPON( XX (  7 ) ,  3 ) ; 
QUEUE( 2) , 1 ; 

ACTIVITY/2, ATRI B( 2) ; 

EVENT, 4; 

ACTIVITY,,, QUE4; 

QUE3  ASS IGN , ATRI 8  <  2 ) -EXPON ( XX! 8 ) , 4 ) ; 
EVENT, 5; 

QUEUE! 3) , 1 ; 

ACTIVITY/3, ATRIB! 2) ; 

EVENT, 6; 

QUE4  ASSIGN, ATRI B ( 2 ) =USERF ( 1 ) ; 

EVENT, 7; 

QUEUE! 4) ,1 ; 

ACTIVITY/4, ATRIB! 2); 

EVENT, 8; 

GOON , 1 ; 

ACTIVITY, ,XX(21 ) ,LOOP; 
ACTIVITY; 

EVENT, 10; 

TERMINATE; 

LOOP  EVENT, 9; 

ACTIVITY,,, QUE1; 

ENDNETWORK; 

INIT, 0,1100; 


COUNT  EXTERNAL  ARRIVALS 

COUNT  ENTITIES  ENTERING  QUEUE  1 

SERVER  #1 
WORK  VARIABLE  #1 

SERVER  #2 
WORK  VARIABLE  #2 

COUNT  ENTITIES  ENTERING  NODE  3 

SERVER  #3 
WORK  VARIABLE  #3 

COUNT  ENTITIES  ENTERING  NODE  4 

SERVER  #4 
WORK  VARIABLE  #4 

COLLECT  SOJOURN  TIME 
COUNT  ENTITIES  LOOPING  BACK 
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FORTRAN  Subroutines  for  Simulation 
XX()  VARIABLE  DEFINITIONS 

XX(1)  x  TRAFFIC  INTENSITY  AT  EACH  NODE  IN  NETWORK 
XX ( 2 )  =  COEF  OF  VAR  OF  SERVICE-TINES  AT  4TH  NODE 

XX(3)  =  1ST  MOMENT  OF  SERVICE-TIMES  AT  4TH  NODE 
XX(4)  x  2ND  MOMENT  OF  SERVICE-TIMES  AT  4TH  NODE 
XX(5)  >  STANDARD  DEV  OF  SERVICE-TIMES  AT  4TH  NODE 

XX(6)  x  MEAN  OF  SERVICE  TIMES  AT  NODE  1 

XX( 7 )  x  MEAN  OF  SERVICE  TIMES  AT  NODE  2 

XX(8 )  »  MEAN  OF  SERVICE  TIMES  AT  NODE  3 

XX( 9 )  S  TIME  TO  BEGIN  COLLECTING  STATISTICS 

XX<10)  x  LAMBDAO  OF  MAX  ENTROPY  DISTRIBUTION 
XX(ll)  >  LAMBDA 1  OF  MAX  ENTROPY  DISTRIBUTION 
XX(12)  =  LAMBDA2  OF  MAX  ENTROPY  DISTRIBUTION 
XX( 13)  x  COEFFICIENT  'C‘  IN  ACCEPTANCE-REJECTION  METHOD 
XX( 14)  x  MEAN  OF  MAJORIZING  DISTRIBUTION 
XX(1S)  x  VARIANCE  OF  MAJORIZING  DISTRIBUTION 

XX( 16)  x  COEFFICIENT  ‘Cl*  OF  HYPEREXPONENTIAL  DISTRIBUTION 
XX( 17)  *  MEAN  OF  FIRST  EXPONENTIAL  DISTRIBUTION 
XX( 18 )  »  COEFFICIENT  *C2*  OF  HYPEREXPONENTIAL  DISTRIBUTION 
XX( 19)  >  MEAN  OF  SECOND  EXPONENTIAL  DISTRIBUTION 

XX(20)  x  PROBABILITY  OF  AN  ENTITY  COING  FROM  NODE  1  TO  NODE  3 
XX( 21 )  x  PROBABILITY  OF  AN  ENTITY  GOING  FROM  NODE  4  TO  NODE  1 
XX (22)  x  TWICE  MEAN  NUMBER  IN  QUEUE  «4  (FOR  QUANTILE  ESTIMATION) 

XX (23)  x  NUMBER  OF  ENTITIES  COMPLETING  NETWORK 
XX( 24)  x  SUM  OF  SOJOURN  TIMES 
XX (25)  x  AVERAGE  SOJOURN  TIME 

XX( 26)  x  FLAG  INDICATING  NNQ(4)  >  XX(22)  ( 0=FALSE, 1=TRUE) 

XX( 27 )  >  LAST  TIME  WHEN  XX(26)  WAS  SET  TO  1 
XX( 28 )  x  SUM  OF  TIME  INTERVALS  WHEN  XX(26)»1 
XX( 29 )  x  PROPORTION  OF  TIME  WHEN  XX(26)xl 

XX( 30 )  x  NUMBER  OF  EXTERNAL  ARRIVALS  TO  NODE  1 
XX (31)  x  EXTERNAL  ARRIVAL  RATE  TO  NODE  1 

XX( 32)  x  NUMBER  OF  ARRIVALS  TO  NODE  1 

XX(  33)  x  NUMBER  OF  ARRIVALS  TO  NODE  3 

XX( 34)  x  STANDARDIZED  ROUTING  VARIABLE  (1,3) 

XX( 35 )  -  NUMBER  OF  ARRIVALS  TO  NODE  4 

XX ( 36)  >  NUMBER  OF  ARRIVALS  TO  NODE  1  FROM  NODE  4 

XX( 37)  x  STANDARDIZED  ROUTINC  VARIABLE  (4,1) 
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XX (39)  *  TOTAL  NUMBER  OF  SERVICE  COMPLETIONS  AT  ALL  NODES 

XX(40)  =  NUMBER  OF  SERVICE  COMPLETIONS  AT  NODE  1 

XX(41 )  *  INTERMEDIATE  SUM  FOR  WORK  VARIABLE  #1 

XX (42)  »  STANDARDIZED  WORK  VARIABLE  #1 

XX (43)  -  NUMBER  OF  SERVICE  COMPLETIONS  AT  NODE  2 

XX (44)  =  INTERMEDIATE  SUM  FOR  WORK  VARIABLE  #2 

XX( 45 )  *  STANDARDIZED  WORK  VARIABLE  #2 

XX(46)  =  NUMBER  OF  SERVICE  COMPLETIONS  AT  NODE  3 

XX( 47 )  ®  INTERMEDIATE  SUM  FOR  WORK  VARIABLE  #3 

XX(48 )  -  STANDARDIZED  WORK  VARIABLE  #3 

XX(49)  *  NUMBER  OF  SERVICE  COMPLETIONS  AT  NODE  4 

XX(SO)  *  INTERMEDIATE  SUM  FOR  WORK  VARIABLE  #4 

XX (51)  *  STANDARDIZED  WORK  VARIABLE  #4 

XX(52)  =  SUM  OF  SERVICE  TIMES  AT  NODE  *1 

XX (53)  =  SUM  OF  SQUARED  SERVICE  TIMES  AT  NODE  #1 

XX(54)  =  SQUARED  COEF  OF  VAR  OF  SERVICE  TIMES  AT  NODE  #1 

XX(SS)  3  SUM  OF  SERVICE  TIMES  AT  NODE  #2 

XX (56)  3  SUM  OF  SQUARED  SERVICE  TIMES  AT  NODE  #2 

XX ( S 7 )  3  SQUARED  COEF  OF  VAR  OF  SERVICE  TIMES  AT  NODE  #2 

XX(58 )  -  SUM  OF  SERVICE  TIMES  AT  NODE  #3 

XX(59 )  3  SUM  OF  SQUARED  SERVICE  TIMES  AT  NODE  #3 

XX(60)  3  SQUARED  COEF  OF  VAR  OF  SERVICE  TIMES  AT  NODE  #3 

XX (61)  3  SUM  OF  SERVICE  TIMES  AT  NODE  #4 

XX (62)  3  SUM  OF  SQUARED  SERVICE  TIMES  AT  NODE  #4 

XX( 63)  3  SQUARED  COEF  OF  VAR  OF  SERVICE  TIMES  AT  NODE  #4 

PROGRAM  MAIN 
DIMENSION  NSET(IOOOO) 

COMMON/SCOMl/ATRIB( 100) ,DD( 100) ,DDL( 100 ) , DTNOW, II ,MFA,MSTOP,NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE , SS ( 1 00 ) , SSL( 1 00 ) , TNEXT , TNOW , XX ( 1 00 ) 
COMMON  QSET( 10000) 

EQU I VALENCE ( NSET( 1 ) ,QSET( 1 ) ) 

NNSET= 10000 
NCRDR=5 
NPRNT=6 
NTAPE® 7 
N PLOT 3 2 

OPEN ( UN I T= 1 , F I LE= ’ QN A . I N ' , S TATUS  * ' NEW  * ) 

OPEN(UNIT32,FILE=’ JACKSON. IN* .STATUS® ' NEW* ) 
OPEN(UNIT=3,FILE3'RESPONSE.OUT’ ,STATUS=’NEW’ ) 
OPEN(UNIT34,FILE-,ROUTINC.OUT' .STATUSs’NEW* ) 

OPEN (UNIT=8,FILE=' WORK. OUT’ ,STATUS=,NEW’ ) 

II  ®  1000 
WRITE  (1,1)  II 
WRITE  (2,1)  II 
CALL  SLAM 
STOP 

1  FORMAT  ( 3X,I6) 

END 


A. 3 


n  nn  o  r>  o  non  o  o  r> 


SUBROUTINE  INTLC 

COMMON /S COM 1 /ATRI B ( 100) ,DD( 100) ,DDL(100) ,DTN0W, II ,MFA,MSTOP,NCLNR 
1 ,  NCRDR ,  NPRNT ,  NNRUN ,  NNSET ,  HTAPE  , SS (1 00 ) , SSL(1 00 ) ,  TNEXT ,  TNOW ,  XX  (1 00 ) 
XX ( 5 )  =  XX<4)  -  XX(3)*XX(3) 

IF  (XX(2) .GT. 1 . )  THEN 

XX( 16)  *  0.5  +  SQRT(0.25  -  0.5*XX( 3 ) *XX ( 3 ) /XX ( 4 ) ) 

XX( 17)  *  0.5*XX( 3)/XX(16 ) 

XX (18)  *  1.0  -  XX( 16) 

XX( 19 )  »  0 .5*XX( 3)/XX( 18 ) 

END  IF 

DO  10  1*23 ,  63 
XXd)  =  0.0 
10  CONTINUE 
RETURN 
END 


FUNCTION  USERFd) 

COMMON/SCOM1 /ATRI B  d  00 ) , DD(1 00 ) , DDL(1 00 ) , DTNOW , 1 1 , MFA , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN ,  NNSET ,  NTAPE ,  SS  (1 00 ) ,  SSLd 00 ) , TNEXT , TNOW , XX (1 00 ) 
INTEGER  ACCEPT 
REAL  F,  G,  PI,  SDEV,  U,  Y 
PARAMETER  (PI  *  3.1415927) 

SELECT  DISTRIBUTION  FOR  SERVICE  TIMES 

IF  (XX(2).LT.l.)  THEN 

GENERATE  MAX  ENTROPY  DISTRIBUTION  USING  ACCEPT-REJECT 

ACCEPT  =  0 

SDEV  *  SORT ( XX ( 15) ) 

10  IF  (ACCEPT. EQ.O)  THEN 

Y  *  RNORM(XX( 14) , SDEV, 9) 

IF  (Y.GE.O.)  THEN 
U  *  ORANDdO) 

F  *  EXP(-1.  -  XX(10)  -  XXd  1  )*Y  -  XX(12)*Y**2) 

G  *  (XX(  13)/SQRT( 2*XX( 15 )*PI ) )* 

&  EXP( -( 1/(2. *XX ( 15) ) )*( Y-XX( 14) )**2) 

IF  (U.LE.F/G)  ACCEPTS 
END  IF 
GO  TO  10 
END  IF 
USERF  =  Y 
ELSE 

GENERATE  HYPEREXPONENTIAL  DISTRIBUTION  USING  COMPOSITION 

U  =  DRAND( 9 ) 

IF  (U.LE.XX( 16) )  THEN 

USERF  *  EXPON(XX( 17 ) , 10) 

ELSE 
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USERF  =  EXPON(XX( 19 ) , 10 ) 
END  IF 
END  IF 
RETURN 
END 


SUBROUTINE  EVENT ( I ) 

C0MM0N/SC0H1/ATR IB(100),DD(100),DDL(100), DTNOW ,11, MFA , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE , SS ( 100 ) , SSL( 100 ) , TNEXT , TNOW , XX ( 1 00 ) 
IF  (TNOW.LE.XX< 9 ) )  RETURN 
CO  TO  (1,  2,  3,  4,  5,  6,  7,  8,  9,  10),  I 

COUNT  EXTERNAL  ARRIVALS  TO  NODE  1 

1  XX ( 30 )  =  XX( 30 )  +  1. 

RETURN 

COUNT  ENTITIES  ARRIVING  TO  NODE  1 

2  XX( 32)  =  XX( 32)  +  1. 

RETURN 

COLLECT  WORK  VARIABLE  #1  DATA 

3  XX (40)  *  XX(40)  +  1. 

XX ( 4 1 )  »  (ATRIB(2)  -  XX(6))/XX(6) 

XX(52)  =  XX(52 )  +  ATRIBC2) 

XX ( S3)  =  XX(53 )  ♦  ATRIB( 2 )*ATRIB( 2 ) 

RETURN 

COLLECT  WORE  VARIABLE  #2  DATA 

4  XX(43)  =  XX(43)  +  1. 

XX (44)  *  (ATRIB(2)  -  XX ( 7 )  )/XX( 7 ) 

XX (55)  =  XX(55)  +  ATRIB(2) 

XX (56)  =  XX ( 5  6 )  +  ATRIB(2)*ATRIB(2) 

RETURN 

COUNT  ENTITIES  ARRIVING  TO  NODE  3 

5  XX ( 33)  =  XX( 33)  +  1. 

RETURN 

COLLECT  WORK  VARIABLE  #3  DATA 

6  XX (46)  *  XX( 46 )  +  1. 

XX(47)  »  ( ATRIB(2)  -  XX(8))/XX(8) 

XX (58)  *  XX(58)  +  ATRIB(2) 

XX(59)  *  XX(S9)  +  ATRIB( 2)*ATRIB( 2 ) 

RETURN 
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COUNT  ENTITIES  ARRIVINC  TO  NODE  4,  AND 
SET  FUG  IF  NNQI4)  >  XXC2),  STORE  TNOW 

7  XXI 35)  =  XXI 35)  +  1. 

IF  I INNQI4 )+l .GT.XXI22) ) .AND.  (XXI26 )  .EQ.  0 . ) )  THEN 
XXI 26)  *  1. 

XX( 27)  =  TNOW 
END  IF 
RETURN 

COLLECT  WORK  VARIABLE  #4  DATA,  AND 
COLLECT  TIME  INTERVAL  NNQI4)  >  XXI 22) 

8  XXI49)  =  XXI 49 )  +  1. 

XX (50)  =  (ATRIBI2)  -  XX(3))/XX(5) 

XXI 61 )  =*  XXI 61 )  ♦  ATRIBI2) 

XXI 62 )  *  XXI 62 )  *  ATRI B( 2 ) *ATRI B( 2 ) 

IF  ( (NNQ(4) . LT.XXI 22) ) .AND. ( XX ( 26 ) . EQ . 1 . ) )  THEN 
XXI 26)  =  0. 

XXI 28)  =  XXI 28)  +  TNOW  -  XXI27) 

END  IF 
RETURN 

COUNT  ENTITIES  LOOPINC  BACK  TO  NODE  1 

9  XXI 36)  *  XXI 36)  +  1. 

RETURN 

COLLECT  DATA  FOR  AVERAGE  SOJOURN  TIME 

10  XXI23)  =  XXI 23)  +  1. 

XXI 24)  *  XXI 24)  +  TNOW  -  ATRIBI 1 ) 

RETURN 

END 


SUBROUTINE  OTPUT 

COMMON/SCOM 1 /ATR IB(100),DD(100),DDL(100)r  DTNOW , 1 1 , MF A , MSTOP , NCLNR 
1 , NCRDR , NPRNT , NNRUN , NNSET , NTAPE , SS ( 100 ) , SSL! 1 00 ) , TNEXT , TNOW , XXI 1 00 ) 
INTEGER  METHOD,  NNODES,  OFTIONI5),  TYPE 
REAL  RATEI4),  ROUTE! 4, 4),  SERVICE! 4) 

CALCULATE  AVERAGE  SOJOURN  TIME 

XXI 25 )  *  XX(24)/XX(23) 

CALCULATE  PROPORTION  OF  TIME  NNQI4)  EXCEEDED 
TWICE  THE  MEAN  NUMBER  IN  QUEUE  FROM  QNA 

IF  ( (NNQI4 ) .LT.XXI 22 ) ) .AND. (XX( 26 ) . EQ. 1 . ) )  THEN 
XXI26)  *  0. 

XXI 28)  *  XXI 28)  +  TNOW  -  XX! 27) 
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END  IF 

XX(29)  =  XX(28)/1000 

CALCULATE  OBSERVED  EXTERNAL  ARRIVAL  RATE 
XX( 31 )  =  XX( 30)/1000 

CALCULATE  STANDARDIZED  ROUTING  CONTROL  VARIABLES 

XX( 34 )  =  ( XX (33)  -  XX ( 32)*XX( 20) )/SQRT(XX( 32)*( 1 . -XX(20) )*XX( 20 ) ) 
XX( 37 )  =  ( XX< 36 )  -  XX ( 35 )*XX( 21 ) )/SQRT(XX( 35 )*( 1 . —XX ( 2 1 ) )*XX( 21 ) ) 

CALCULATE  STANDARDIZED  WORE  VARIABLES 

XX( 39)  =  XX( 40 )  +  XX(43)  +  XX(46)  +  XX(49) 

XX ( 42)  =  3*SQRT( XX(40 ) )/XX( 39 )*XX(4i ) 

XX(45)  *  3*SQRT(XX( 43) )/(XX( 39 )*! 1 .-XX( 33 ) /XX ( 32) ) )*XX( 44 ) 

XX(48)  =  3*SQRT( XX! 46 ) )/( XX( 39 )*XX( 33)/XX( 32) )*XX( 47 ) 

XX (51)  =  3*SQRT(XX( 49 ) )/XX( 39 )*XX(41 ) 

CALCULATE  OBSERVED  MEAN  SERVICE  TIMES 

SERVICE(l)  =  XX(52)/XX(40) 

SERVICE! 2)  =  XX( 55 )/XX( 43 ) 

SERVICE! 3)  =  XX(5B)/XX(46) 

SERVICE(4)  =  XX(6I)/XX!49) 

CALCULATE  SQUARED  COEFFICIENT  OF  VARIATION  OF  SERVICE  TIMES 

XX(54)  =»  ( (XX(40)**2)*XX(53)  -  XX(40)*(XX(52)**2) )/ 
&(!XX(40)-1)*XX(52)**2) 

XXI57)  =  (!XX(43)**2)*XX!56)  -  XX!43)*(XX(55 )**2) )/ 

&( (XX(43) -I )*( XX!55 )**2) ) 

XX( 60)  =  ! (XX(46 )**2)*XX!59 )  -  XX!46 )*(XX(58 )**2 ) )/ 

&( (XX(46)-1)*(XX(58)**2)) 

XX! 63)  =  ! (XX! 49 )**2)*XX! 62)  -  XX(49 )*!XX( 61 )**2) )/ 

&! (XX(49 )-l )*<XX( 61 )**2 ) ) 

WRITE  TO  'QNA. IN' 

METHOD  *  3 
NNODES  =  4 
OPTION! 1)  *  1 
OPTION! 2)  *  2 
OPTION! 3)  =  0 
OPTION! 4)  =  -1 
OPTION! 5)  =  1 
TYPE  *  1 
ROUTE! 1,1)  =  0. 

ROUTE! 1,2)  =  1.  -  XX! 33)/XX( 32) 

ROUTE! 1,3)  =  XX(33)/XX( 32) 

ROUTE! 1,4)  *  0. 

ROUTE! 2,1)  =  0. 
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ROUTE(2,2)  =  0. 

R0UTE(2, 3)  =  0. 

ROUTE (2, 4)  =  1. 

R0UTE( 3,1)  =  0. 

R0UTE( 3,2)  =  0. 

R0UTE( 3, 3)  =  0. 

R0UTE( 3,4)  =  1. 

R0UTE(4,1)  =  XX ( 36) /XX ( 35 ) 

ROUTE(4,2)  =  0. 

R0UTE(4, 3)  =  0. 

ROUTE(4,4)  =  0. 

RATE(l)  =  XX( 31 ) 

RATE(2)  =  0. 

RATE( 3)  =  0. 

RATE(4)  =  0. 

WRITE  (1,1)  METHOD 
WRITE  (1,2)  NNODES,  TYPE 

WRITE  (1,3)  OPTION(l),  0PTI0N(2),  0PTI0N(3),  0PTI0N(4),  0PTI0N(5) 
DO  10  1=1,4 

WRITE  (1,4)  R0UTE( I ,1 ) ,  R0UTE(I,2),  R0UTE(I,3),  R0UTE(I,4) 

10  CONTINUE 

WRITE  (1,4)  RATE ( 1 ) ,  RATE( 2) ,  RATE( 3) ,  RATE(4) 

WRITE  (1,4)  SERVICE(l),  SERVICE(2),  SERVICE(3),  SERVICED) 

WRITE  (1,4)  XX(54),  XX(57),  XX(60),  XX(63) 

WRITE  TO  ‘JACKSON. IN' 

OPTION(l)  =  1 
OPTION(2)  =  0 
OPTION( 3)  =  0 
OPTION(4)  =  -1 
OPTION(5 )  =  1 
WRITE  (2,1)  METHOD 
WRITE  (2,2)  NNODES ,  TYPE 

WRITE  (2,3)  OPTION(l),  OPTION(2),  OPTION(3),  OPTION(4),  OPTION(5) 
DO  20  1=1,4 

WRITE  (2,4)  ROUTE ( 1,1),  R0UTE(I,2),  R0UTE(I,3),  ROUTE(I,4) 

20  CONTINUE 

WRITE  (2,4)  RATE(l),  RATE (2 ) ,  RATE( 3) ,  RATE(4 ) 

WRITE  (2,4)  SERVICE(l),  SERVICED),  SERVICED),  SERVICE(4) 

WRITE  TO  ‘RES PONSE.OUT* 

NUMBER  OF  OBSERVATIONS,  AVE  SOJOURN  TIME,  AND  P(NNQ(4)>2*EN) 

WRITE  (3,5)  XX{23),  XX(25),  XX(29) 

WRITE  TO  ‘ROUTING. OUT’ 

ROUTING  CONTROLS  R13,  R14 

WITH  RESPECTIVE  NUMBER  OF  OBSERVATIONS  ON  EACH 

WRITE  (4,4)  XX( 32 ) ,  XX(34),  XX(35),  XX(37) 

C 

C  WRITE  TO  ‘WORK. OUT’ 
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WORK  VARIABLES  Wl,  W2,  W3,  W4 
WITH  HUMBER  OF  OBSERVATIONS 

WRITE  (8,6)  XX( 39) ,  XX(42),  XX(45),  XX(48),  XX(51) 
RETURN 

1  FORMAT(2X,I2) 

2  FORMAT ( 2 ( 2X , I 2 ) ) 

3  FORMAT! 5 (2X, 12)) 

4  FORMAT (4(2X,E13.6)) 

5  FORMAT! 3(2X,E13.6) ) 

6  FORMAT ( 5 ( 2X , El  3 . 6 ) ) 
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FORTRAN  Code  for  Proeraa  CONTROL 
PROGRAM  CONTROL 

THIS  PROGRAM  CALCULATES  THE  VARIANCE  RATIO  (VR)  AND 
THE  95X  CONFIDENCE  LIMITS  (LCL)  AND  (UCL)  ON  THE 
CONTROLLED  RESPONSES. 

THE  INPUT  DATA  IS  SUPPLIED  BY  SEVEN  FILES: 

(1)  RESPONSE. OUT  CONTAINS: 

NUMBER  OF  OBSERVATIONS  OF  SOJOURN  TIME, 
SOJOURN  TIME  (UNCONTROLLED  RESPONSE),  AND 
P(NNQ(4)>2*EN)  (UNCONTROLLED  RESPONSE). 

(2)  ROUTING. OUT  CONTAINS: 

NUMBER  OF  OBSERVATIONS  ON  R13, 

R13  (STANDARDIZED  ROUTING  CONTROL), 

NUMBER  OF  OBSERVATIONS  ON  R41,  AND 
R41  (STANDARDIZED  ROUTING  CONTROL). 

(3)  WORK. OUT  CONTAINS: 

TOTAL  NUMBER  OF  SERVICE  COMPLETIONS,  AND 
Wl,  W2,  W3,  W4  (STANDARDIZED  WORK  VARIABLES). 

(4)  JACKSON. OUT  CONTAINS  EXTERNAL  CONTROL  VARIATES 
USING  SHARON'S  (1986)  METHOD  AND  M/M/1  FORMULAS: 

SOJOURN  TIME,  AND 
WAITING  TIME  IN  4TH  QUEUE. 

(5)  QNA.OUT  CONTAINS  EXTERNAL  CONTROL  VARIATES 
USING  SHARON'S  (1986)  METHOD  AND  G/G/l  FORMULAS: 

SOJOURN  TIME,  AND 
WAITING  TIME  IN  4TH  QUEUE. 

(6)  Ml 6 .OUT  CONTAINS  THE  MEANS  OF  THE  EXTERNAL 
CONTROL  VARIATES  USING  M/M/1  FORMULAS. 

(7)  G16.0UT  CONTAINS  THE  MEANS  OF  THE  EXTERNAL 
CONTROL  VARIATES  USING  G/G/l  FORMULAS. 

THE  OUTPUT  DATA  IS  WRITTEN  TO  FOUR  FILES  ( SOJOURN. VR, 
SOUJOURN.CI ,  QUANTILE. VR  AND  QUANTILE. Cl)  IN  THE 
FOLLOWING  FORMAT: 

THE  *.VR  FILES  CONTAIN: 

MEAN  OF  UNCONTROLLED  RESPONSE 
MEAN  OF  CONTROLLED  RESPONSE 
VARIANCE  RATIO 
THE  *.CI  FILES  CONTAIN: 

9SX  LOWER  CONFIDENCE  LIMIT 
95 X  UPPER  CONFIDENCE  LIMIT 

EACH  RECORD  IN  THE  OUTPUT  FILES  ARE  THE  RESULT  OF 
TWENTY  INPUT  RECORDS  (I.B.  TWENTY  SIMULATION  RUNS). 

EVERY  TENTH  RECORD  IN  THE  OUTPUT  FILES  REPRESENTS  A 
REPLICATION  OF  THE  RESULTS  WITH  ONE  CONTROL  VARIATE. 

INPUT  VARIABLES 

Y( I ,K)  =  TWO  RESPONSE  VARIABLES  (20  OB S.  EA.) 

C( J,K)  -  TEN  CONTROL  VARIABLES  (20  OBS.  EA.) 
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CC1(2,16)  a  MEANS  OF  ANALYTIC  CONTROLS  USING  QNA  (C/G/l) 

MM1(2, 16)  =  MEANS  OF  ANALYTIC  CONTROLS  USING  QNA  (M/M/1) 

OBSY(K)  =  NUMBER  OF  OBSERVATIONS  IN  SIMULATION  MODEL 

OF  THE  SOJOURN  TIME  TO  YIELD  *KTH*  OBSERVATION 
OF  AVERAGE  SOJOURN  TIME 

OBSR( J,K)  »  NUMBER  OF  OBSERVATIONS  IN  SIMULATION  MODEL 

OF  THE  ROUTE  TAKEN  TO  COMPUTE  THE  'KTH'  OBSERVATION 
OF  THE  ROUTING  CONTROL  VARIABLES  R13  AND  R41 
OBSW(K)  =  TOTAL  NUMBER  OF  SERVICE  COMPLETIONS  FOR  THE 
•KTH*  OBSERVATION  OF  THE  FOUR  WORK  VARIABLES 

REAL  Y(2,20),  C(10,20),  GG1(2,16),  MM1(2,16) 

REAL  OBSY(20) ,  OBSR(2,20),  OBSW(20) 

OUTPUT  VARIABLES 

YBAR(I)  a  GRAND  UNCONTROLLED  MEAN  OF  THE  TWO  RESPONSES 
YCBAR( I , J)  >  CRAND  CONTROLLED  MEAN  OF  THE  TWO  RESPONSES 
USING  THE  ' JTH'  CONTROL 

VR( I , J)  a  VARIANCE  RATIO  OF  THE  ' ITH'  RESPONSE  US  INC  THE 
•JTH*  CONTROL 

LCL( I , J)  a  95X  LOWER  CONFIDENCE  LIMIT  FOR  THE  CONTROLLED 

MEAN  OF  THE  'ITH*  RESPONSE  USING  THE  'JTH'  CONTROL 
UCL( I , J)  a  9SX  UPPER  CONFIDENCE  LIMIT  ... 

REAL  YBAR(2),  YCBAR(2,10),  VR(2,10),  LCL(2,10),  UCL(2,10) 

INTERMEDIATE  VARIABLES  USED  IN  CALCULATIONS 

BHAT ( I ,  J )  a  ESTIMATED  OPTIMAL  CONTROL  COEFFICIENT 

FOR  THE  'ITH'  RESPONSE  AND  'JTH'  CONTROL 
CBAR(J)  =  MEAN  OF  THE  'JTH'  CONTROL  VARIABLE 
COVYC( I ,J)  =  COVARIANCE  BETWEEN  'ITH'  RESPONSE  AND 
'JTH'  CONTROL 

D2(J)  =  D  SQUARED  STATISTIC  OF  'JTH'  CONTROL 

I , J,K  a  ITERATION  VARIABLES/ INDEX  VARIABLES 

MU( J)  a  KNOWN  MEAN  OF  'JTH'  CONTROL 

POINT  a  EXPERIMENTAL  DESIGN  POINT  (INDEX  VARIABLE) 

REP  a  ITERATION  VARIABLE  (REPLICATION) 

REPS  a  TOTAL  NUMBER  OF  REPLICATIONS  OF  EXPERIMENT,  WHERE 
AN  EXPERIMENT  IS  20  OBSERVATIONS  OF  THE  SIMULATION 
MODEL  UNDER  THE  SAME  INTIAL  CONDITIONS 
RHOYC( I , J)  =  PEARSON'S  PRODUCT-MOMENT  CORRELATION  STATISTIC 
BETWEEN  'ITH'  RESPONSE  AND  'JTH'  CONTROL 
S2C( J)  a  SAMPLE  VARIANCE  OF  'JTH*  CONTROL 
S2Y( I )  «  SAMPLE  VARIANCE  OF  'ITH*  RESPONSE 
S2YC(I,J)  a  S  SQUARED  (Y  DOT  C)  STATISTIC  OF  'ITH' 

RESPONSE  AND  'JTH'  CONTROL 
SUMC(J)  »  SUM  OF  OBSERVATIONS  ON  'JTH*  CONTROL 
SUMC2(J)  a  SUM  OF  SQUARED  OBSERVATIONS  ON  ’JTH’  CONTROL 
SUMY ( I )  a  SUM  OF  OBSERVATIONS  ON  'ITH’  RESPONSE 
SUMY2(I)  >  SUM  OF  SQUARED  OBSERVATIONS  ON  'ITH'  RESPONSE 
TOTOBS  >  TOTAL  NUMBER  OF  OBSERVATIONS 
TOTOBSR(I)  «  TOTAL  NUMBER  OF  SIMULATION  OBSERVATIONS  TO 

PRODUCE  THE  20  OBSERVATIONS  OF  THE  TWO  ROUTING 
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C  CONTROLS 

C  TOTOBSW  *  TOTAL  NUMBER  OF  SIMULATION  OBSERVATIONS  OF 

C  SERVICE  COMPLETIONS  USED  IN  COMPUTING  THE  20 

C  OBSERVATIONS  ON  THE  FOUR  WORK  VARIABLES 

C  TOTOBSY  *  TOTAL  NUMBER  OF  SIMULATION  OBSERVATIONS  TO  PRODUCE 

C  THE  20  OBSERVATIONS  OF  AVERAGE  SOJOURN  TIME 

C  WC( J,K)  *  WEIGHT  OF  ’KTH'  OBSERVATION  ON  ' JTH’  CONTROL 

C  WY(I,K>  =  WEIGHT  OF  ’KTH’  OBSERVATION  ON  ' ITH*  RESPONSE 

C  ZERO  =  INITIALIZATION  VARIABLE  <»  0.0) 

C 

REAL  BHAT(2,10),  CBAR(IO),  C0VYC(2,10),  D2(10) 

REAL  MU( 10) ,  RHOYC(2,10),  S2C(10),  S2Y(2),  S2YC(2,10) 

REAL  SUMC(IO),  SUMC2U0),  SUMY(2)f  SUMY2(2) 

REAL  TOTOBSR(2),  TOTOBSY,  TOTOBSW 
REAL  WC(10,20),  WY{2,20),  ZERO 
INTEGER  I,  J,  K,  POINT,  REP,  REPS,  TOTOBS 
C 

C  INTERACTIVE  USER  INPUT 

C 

WRITE  (5,*)  'ENTER  TOTAL  NUMBER  OF  OBSERVATIONS ’ 

READ  (5,*)  TOTOBS 
REPS  *  TOTOBS/20 

WRITE  (5,*)  ’ENTER  THE  DESIGN  POINT  NUMBER  (1  TO  16)' 

READ  (5,*)  POINT 
C 

C  READ  ANALYTIC  RESULTS 

C 

OPEN  (UNIT-1, FILE* 'C16. OUT’, STATUS* 'OLD') 

OPEN  (UNIT=2, FILE* ’M16. OUT* , STATUS* ’OLD’ ) 

DO  10  1-1,16 

READ  (1,1)  CC1(1,I),  GG1 (2,1) 

READ  (2,1)  MMl(l,I),  MM1(2,I) 

10  CONTINUE 

1  FORMAT  (2(2X,E13.6) ) 

CLOSE  (1) 

CLOSE  (2) 

C 

C  INITIALIZE  KNOWN  MEANS  OF  CONTROLS 
C 

ZERO  *  Q.0E+00 
DO  20  1*1,6 
MU( I )  =  ZERO 
20  CONTINUE 

MU<7)  *  M!l(l, POINT) 

MU( 8 )  *  MM1(2, POINT) 

MU(9)  -  GG 1(1, POINT) 

MU(10)  *  GC1(2, POINT) 

C 

C  EXPERIMENTAL  DATA  FILES  FOR  INPUT 
C 

OPEN  (UNIT-1 , FILE* ’ RESPONSE. OUT’ , STATUS* 'OLD ’ ) 

OPEN  (UNIT-2, FILE- ’ROUTING. OUT’ , STATUS* ’OLD’ > 

OPEN  (UNIT- 3, FILE* ’WORK. OUT’ .STATUS* ’OLD’ ) 
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OPEH  (UNIT=4,FILE=‘ JACKSON. OUT'  ,STATUS-'OLD' ) 

OPEN  (UNIT-8, FILE- *QNA. OUT* ,STATUS-'OLD' ) 

OUTPUT  DATA  FILES 

OPEN  ( UNIT-9, FILE- ‘ SOJOURN . VR * , STATUS- 'NEW ) 

OPEN  (UNIT-10, FILE- 'SOJOURN. Cl* ,STATUS='NEW ) 

OPEN  (UNIT- 11, FILE-' QUANTILE. VR' , STATUS- 'NEW ) 

OPEN  (UNIT- 12, FILE- 'QUANTILE. Cl* .STATUS- 'NEW' ) 

DO  1000  REP- 1, REPS 

INITIALIZE  SUMMATION  VARIABLES 

TOTOBSY  -  ZERO 
TOTOBSW  -  ZERO 
TOTOBSR(l)  -  ZERO 
TOTOBSR( 2)  -  ZERO 

READ  EXPERIMENTAL  DATA 

DO  30  J-1,20 

READ  (1,2)  OBSY(J),  Y(1,J),  Y(2,J) 

TOTOBSY  -  TOTOBSY  +  OBSY(J) 

READ  (2,3)  OBSR(l.J),  C(1,J),  OBSR(2,J),  C(2,J) 
TOTOBSR(l)  -  TOTOBSR(l)  *  OBSR(l,J) 

TOTOBSR( 2)  »  TOTOBSR(2)  +  OBSR(2,J) 

READ  (3,4)  OBSW(J),  (C(I,J),  1=3,6) 

TOTOBSW  -  TOTOBSW  +  OBSW(J) 

READ  (4,5)  C( 7 , J) ,  C(8,J) 

READ  (8,5)  C(9, J) ,  C(10,J) 

30  CONTINUE 

2  FORMAT  ( 3(2X,E13. 6) ) 

3  FORMAT  (4(2X,E13.6) ) 

4  FORMAT  (5(2X,E13.6) ) 

5  FORMAT  (2( 2X.E13.6) ) 

INITIALIZE  VARIABLES  FOR  SAMPLE  MEANS, 

VARIANCES,  AMD  COVARIANCES 

DO  31  1-1,2 

YBAR(I)  -  ZERO 
SUMY ( I )  =  ZERO 
SUMY2(I)  =  ZERO 
S2Y(I)  =  ZERO 

31  CONTINUE 

DO  32  1-1,10 

CBAR(I)  -  ZERO 
SUMC(I)  -  ZERO 
SUMC2( I )  -  ZERO 
S2C( I )  >  ZERO 
COVYC(l.I)  =  ZERO 
COVYC(2, I )  -  ZERO 


A. 13 


nr>o  o  n  o  non  non 


32 


CONTI HUG 


CALCULATE  WEIGHTS  OF  OBSERVATIONS 
DO  35  J*l,20 

WY ( 1 , J )  =  20 . 0*OBSY ( J )/T0T0BSY 
WY ( 2 , J )  *  1.0 

WC(1,J)  *  20 . 0*OBSR ( 1 , J ) /TOTOBSR ( 1 ) 

WC( 2, J)  =  20 . 0*OBSR( 2 , J ) /T0T0BSR( 2 ) 

DO  33  1*3,6 

WC(I,J)  *  20 . 0*OBSW ( J ) /TOTOBSW 

33  CONTINUE 
DO  34  1=7,10 

WC(I,J)  =  20 . 0*OBS Y ( J ) /TOTOBS Y 

34  CONTINUE 

35  CONTINUE 

CALCULATE  WEIGHTED  SUMS 

DO  42  J=1 ,20 
DO  40  1=1,2 

SUMY ( I )  =  SUMY ( I )  +  WY(I,J)*Y(I,J) 

SUMY2( I )  =  SUMY2( I )  *  WY( I , J)*Y( I , J)**2 

40  CONTINUE 
DO  41  1=1,10 

SUMC(I)  =  SUMC(I)  +  WC(I,J)*C(I,J) 

SUMC2( I )  =  SUMC2( I )  f  WC(I,J)*C(I,J)**2 

41  CONTINUE 

42  CONTINUE 

CALCULATE  WEIGHTED  AVERAGES  AND  WEIGHTED  SAMPLE  VARIANCES 
DO  50  1=1,2 

YBAR(I)  =  SUMY( I )/20 . 0 

52Y( I )  =  (20. 0*SUMY2{ I )  -  SUMY( I )**2)/380 . 0 

50  CONTINUE 
DO  51  1=1,10 

CBAR(I)  =  SUMC( I )/20 .0 

S2C( I)  =  ( 20 . 0*SUMC2( I )  -  SUMC( I )**2)/380.0 

51  CONTINUE 

CALCULATE  SAMPLE  COVARIANCES 

DO  62  1=1,2 
DO  61  J= 1,20 
DO  60  K=1 , 10 

COVYC( I ,K)  =  C0VYC( I ,R>  +  (C(R,J)  -  CBAR(K) )* 
&(Y( I ,J)  -  YBAR( I ) ) / 1 9 

60  CONTINUE 

61  CONTINUE 

62  CONTINUE 
C 

C  ESTIMATE  OPTIMAL  CONTROL  COEFFICIENT  AND 
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C 

C 


PEARSON'S  PRODUCT-MOMENT  CORRELATION  COEFFICIENT 


DO  71  J=l,10 
DO  70  1*1,2 

BHAT( I , J)  =  C0VYC(I,J)/S2C(J) 

RH0YC( I,J)  *  C0VYC{ I , J)/SQRT(S2C( J)*S2Y( I ) ) 

70  CONTINUE 

71  CONTINUE 

CALCULATE  CONTROLLED  MEANS  OF  RESPONSE  VARIABLES 

DO  73  J=l,10 
DO  72  1-1,2 

YCBAR<I,J)  =  YBAR(I)  +  BHAT(I , J)*(CBAR( J)-MU( J) ) 

72  CONTINUE 

73  CONTINUE 

CALCULATE  STATISTICS  'D2'  AND  'S2YC*  TO  ESTIMATE 
VARIANCE  OF  CONTROLLED  MEANS 

VAR[ YCBAR{ I ,  J)  ]  =  D2(I,J)*S2YC(I,J) 

DO  81  J=»l,10 

D2( J)  =  0.05+( 1/19)*(CBAR( J)-MU( J) )**2/S2C( J) 

DO  80  1*1,2 

S2YC( I , J)  =  ( 19/18 )*S2Y( I )*( 1 . C  -  RHOYC( I , J)**2) 

80  CONTINUE 

81  CONTINUE 

ESTIMATE  VARIANCE  RATIOS  AND  95X  CONFIDENCE  LIMITS 

DO  91  J=l,10 
DO  90  1*1,2 

VR( I , J)  =  (D2( J)*S2YC( I , J) )/(S2Y( I )/20 ) 

LCL( I , J)  =  YCBAR( I , J)  -  2 . 101*SQRT(D2( J)*S2YC( I , J) ) 
UCL(I.J)  =  YCBAR(I , J)  +  2 . 101*SQRT(D2( J)*52YC( I , J) ) 

90  CONTINUE 

91  CONTINUE 

OUTPUT  RESULTS 
DO  100  J*l,10 

WRITE  (9,101)  YBAR(l),  YCBAR(1,J),  VR(1,J) 

WRITE  (10,102)  LCL( 1 , J) ,  UCL(1,J) 

WRITE  (11,101)  YBAR(2),  YCBAR(2,J),  VR(2,J) 

WRITE  (12,102)  LCL( 2, J) ,  UCL(2,J> 

100  CONTINUE 

101  FORMAT  (3(2X,E13.6)) 

102  FORMAT  (2(2X,E13.6) ) 

1000  CONTINUE 

END 
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PROGRAM  RESULTS 


THIS  PROGRAM  FINDS  THE  MINIMUM,  MAXIMUM  AND  AVERAGE  OF  THE  THREE 
VARIABLES  (YBAR,  YCBAR,  AND  VR)  FROM  THE  OUTPUT  FILES  OF  THE 
PROGRAM  ’CONTROL1.  IT  ALSO  COMPUTES  THE  COVERAGE  OF  THE  9SX 
CONFIDENCE  INTERVAL  ACAINST  THE  UNCONTROLLED  RESPONSE  (YBAR). 

INPUT  VARIABLES 

Y(1 , J)  *  UNCONTROLLED  MEAN  OF  THE  RESPONSE  VARIABLE 
Y(2,J)  =  CONTROLLED  MEANS  OF  THE  RESPONSE  VARIABLE 
USING  ' JTH'  CONTROL 

Y( 3, J)  =  VARIANCE  RATIO  OBTAINED  USING  ’JTH’  CONTROL 
LCL(J)  *  95X  LOWER  CONFIDENCE  LIMIT  ON  CONTROLLED  MEAN  OF 
THE  RESPONSE  VARIABLE  USING  THE  ’JTH’  CONTROL 
UCL(J)  =  95%  UPPER  CONFIDENCE  LIMIT  ... 

REAL  Y( 3, 10) ,  LCL(IO),  UCL(IO) 

OUTPUT  VARIABLES 

MIN(I,J)  =  MINIMUM  VALUE  OF  *ITH’  INPUT  USING  ’JTH’  CONTROL 
MAX(I,J)  =  MAXIMUM  VALUE  OF  ’ ITH’  INPUT  USING  ’JTH’  CONTROL 
MEAN(I,J)  *  MEAN  VALUE  OF  ’ITH’  INPUT  USING  ’JTH’  CONTROL 
COVER(J)  =  COVERAGE  PROBABILITIES  OF  95X  CONFIDENCE  INTERVAL 
OF  THE  CONTROLLED  RESPONSE  COMPARED  TO  THE 
CRAND  UNCONTROLLED  MEAN 

REAL  MIN( 3,10),  MAX(3,10),  MEAN(3,10),  COVER(IO) 

INTERMEDIATE  VARIABLES  USED  IN  CALCULATIONS 

COUNT ( J )  =  NUMBER  OF  TIMES  THE  MEAN  RESPONSE  FALLS  WITHIN 
THE  CONFIDENCE  INTERVALS 
INFILE  *  INPUT  FILE  TYPE 
I,J  =  ITERATION/INDEX  VARIABLES 
REP  *  ITERATION  VARIABLE  (REPLICATION) 

REPS  =  TOTAL  NUMBER  OF  REPLICATIONS  OF  THE  EXPERIMENT 
L(J)  ■  LABEL  OF  ’JTH’  CONTROL 
N( J)  =  LABEL  OF  ’ITH’  STATISTIC 

INTEGER  COUNT (10),  INFILE,  I,  J,  REP,  REPS 
CHARACTER* 14  L(10) 

CHARACTER* 14  N(3) 

INTIALIZE  MINIMUM  VALUES  TO  A  HIGH  NUMBER 

DATA  (MIN( 1 , J) ,  J*1 , 10)  /10*1 .0E+10/ 

DATA  (MIN( 2, J) ,  J-1,10)  /10*1.0E*10/ 

DATA  (MIN( 3, J) ,  J=l,10)  /10*1.0E+10/ 

INITIALIZE  LABELS 

L(l)  *  *ROUTING( 1,3)  ’ 
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L(2)  *  *ROUTING(  1,4)  * 

L(  3)  =  'WORK( 1 ) 

L(4)  *  ’W0RK(2) 

L(5 )  *  ' WORK( 3) 

L(  6)  =  *W0RR<4) 

L( 7)  =  *  SOJOURN ( M/M/1)' 

L<8>  =  'HAIT4(M/M/1 )  * 

L(9)  =  ' SOJOURN (G/G/l ) * 

L(10)  *  'WAIT4(G/G/1 )  * 

N( 1 )  *  ’YBAR 
N(2)  =  'YBAR(BHAT) 

N( 3)  =  ’VARIANCE  RATIO* 

INTERACTIVE  USER  INPUT 

WRITE  (5,*)  ’ENTER  THE  NUMBER  OF  REPLICATIONS’ 

READ  <5,0  REPS 

WRITE  (5,*)  ’INDICATE  INPUT  FILES  (1=S0J0URN,  2=QUANTILE) ’ 
READ  (S,*)  INFILE 

OPEN  APPROPRIATE  INPUT  DATA  FILES  (CREATED  BY  PROGRAM  ’CVR’) 
IF  ( INFILE. EQ.l)  THEN 

OPEN  (UNIT=1 , FILE* ’SOJOURN. VR’ , STATUS =  ’ OLD ’ ) 

OPEN  ( UNIT=2, FILE* ’SOJOURN. Cl* , STATUS* ’OLD’ ) 

ELSE 

OPEN  ( UNI TO ,FILE=’ QUANTILE . VR ’ , STATUS* ’ OLD ’ ) 

OPEN  ( UNIT*2 , FILE* ’ QUANTILE . Cl 1 , STATUS* ’ OLD ‘ ) 

END  IF 

OPEN  OUTPUT  DATA  FILE 

OPEN  (UNi 1=3, FILE*’ RESULTS. DAT’ ,STATUS=’NEW’ ) 

DO  30  REP* 1, REPS 

READ  INPUT  VARIABLES  FROM  *.VR  FILE 
DO  10  J=l,10 

READ  (1,1)  ( Y( I , J) ,  1*1,3) 

10  CONTINUE 

COMPUTE  MIN,  MAX,  AND  MEAN  OF  THE  INPUT  VARIABLES 

DO  21  1*1,3 
DO  20  J=l,10 

IF  (Y( I , J) .LT.MIN( I, J) )  MI N( I , J)  =  Y(I,J) 

IF  ( Y( I , J) .GT.MAX( I , J) )  MAX ( I , J )  *  Y(I,J) 

MEAN(I,J)  *  MEAN(I,J)  ♦  Y(I,J)/REPS 

20  CONTINUE 

21  CONTINUE 
30  CONTINUE 

C 
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C  READ  CONFIDENCE  LIMITS  FROM  *.CI  FILE 
C 

DO  60  RBP=1 ,REPS 
DO  40  J=1 , 10 

READ  (2,2)  LCL(J),  UCL(J) 

40  CONTINUE 

CALCULATE  COVERAGE  ESTIMATES  OF  CONFIDENCE  INTERVALS 
DO  SO  J=l,10 

IF  ( (LCL( J) . LE.MEAN( 1 , J) ) .AND. (MEAN( 1 , J) . LE.UCL( J) ) )  THEN 
COUNT ( J )  =  COUNT ( J )  +  I 
END  IF 
50  CONTINUE 

60  CONTINUE 

DO  70  J=1 , 10 

COVER ( J )  =  1 . 0*COUNT ( J ) /REPS 
70  CONTINUE 

OUTPUT  RESULTS 

IF  ( INFILE. EQ.l)  THEN 
WRITE  (3,3) 

ELSE 

WRITE  (3,4) 

END  IF 
WRITE  (3,0) 

WRITE  (3,6) 

WRITE  (3,7) 

WRITE  (3,8) 

DO  80  J=l,10 

WRITE  (3,5)  L(J) 

WRITE  (3,9)  N( 1 ) ,  NIN(l.J),  MEAN(1,J),  MAX(1,J),  COVER(J) 

WRITE  (3,11)  N(2) ,  MIN(2,J),  MEAN(2,J),  MAX(2,J) 

WRITE  (3,11)  N( 3) ,  MIN( 3 , J) ,  MEAN(3,J),  MAX(3,J) 

WRITE  (3,8) 

80  CONTINUE 
STOP 

1  FORMAT  ( 3(2X,E13.6) ) 

2  FORMAT  (2(2X,E13.6) ) 

3  FORMAT  ( 15X, 'CONTROL  VARIATE  ANALYSIS  ON  RESPONSE:  SOJOURN  TIME1) 

4  FORMAT  ( 10X, ' CONTROL  VARIATE  ANALYSIS  ON  RESPONSE:  ', 

&  ' PR0B(NNQ(4)  >  2*E[N]) ’ ) 

5  FORMAT  (IX/ US  INC  CONTROL  VARIABLE:  ' ,A14> 

6  FORMAT  ( 3X, '  NAME  MINIMUM  AVERAGE 

&  ’  MAXIMUM  COVERAGE') 

7  FORMAT  ( 3X,  ’ -  -  - -  ', 

4  * . ') 

8  FORMAT  (IX) 

9  FORMAT  ( 3X,A14 , 3(2X,E13. 6) , 3X,F6 . 4) 

11  FORMAT  (3X,A14,3(2X,E13.6)) 

END 
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Appendix  B:  Saaple  Input  Files  to  QNA 


M/M/1  Input  File  at  First 

Design  Point 

1 

(indicates  1  network) 

3 

(indicates  single  customer  class) 

4  1 

(4  nodes,  nonstandard  input) 

100-10 

(input  options) 

0.0  0.6  0.4 

0.0 

(1st  row  of  routing  aatrix) 

0.0  0.0  0.0 

1.0 

(2nd  row  of  routing  aatrix) 

0.0  0.0  0.0 

1.0 

(3rd  row  of  routing  aatrix) 

0.1  0.0  0.0 

0.0 

(4th  row  of  routing  aatrix) 

1.0  0.0  0.0 

0.0 

(external  arrival  rates) 

0.450  0.750 

1.125 

0.450 

(aean  service  tiaes) 

C/C/1  Input  File  at 

First 

Design 

Point 

1 

(indicates  1  network) 

3 

(indicates  single  custoaer  class) 

4  1 

(4  nodes,  nonstandard  input) 

120-10 

(input  options) 

0.0  0.6  0.4 

0.0 

(1st  row  of  routing  aatrix) 

0.0  0.0  0.0 

1.0 

(2nd  row  of  routing  aatrix) 

0.0  0.0  0.0 

1.0 

(3rd  row  of  routing  aatrix) 

0.1  0.0  0.0 

0.0 

(4th  row  of  routing  aatrix) 

1.0  0.0  0.0 

0.0 

(external  arrival  rates) 

0.450  0.750 

1.125 

0.450 

(aean  service  tines) 

1.000  1.000 

1.000 

0.250 

(service  variability  paraaeters) 

B.l 


9 


Appendix  C:  Known  Means  of  External  Control  Variates 


DESIGN  —ANALYTIC  JACKSON  (M/M/1 )-  -ANALYTIC  APPROX.  (G/G/l)  — 


INT 

SOJOURN 

WAIT(4) 

SOJOURN 

WAIT(4) 

i 

0.400000E+01 

0 . 450000E+00 

0. 381177E+01 

0.281197E+00 

2 

0.360000E+02 

0 . 729000E+01 

0 . 329330E+02 

0.455585E+01 

3 

0 . 400000E+01 

0 . 450000E+00 

0.531762E+01 

0.163162E+01 

4 

0 . 360000E+02 

0 . 729000E+01 

0.574688E+-02 

0 . 264290E  +  02 

5 

0.400000E+01 

0.450000E+00 

0. 381173E+01 

0 . 281160E+00 

6 

0 . 360000E+02 

0 . 729000E+01 

0 . 329329E+02 

0.455572E+01 

7 

0 . 400000E+01 

0 . 450000E+00 

0.531792E+01 

0.163188E+01 

8 

0. 360000E+02 

0 . 729000E+01 

0.S74698E+02 

0 . 264300E+02 

9 

0. 417391E+01 

0 . 375000E+00 

0. 398269E+01 

0 . 234173E+00 

10 

0 . 360000E+02 

0 . 607500E+01 

0 . 327806E+02 

0 . 379484E+01 

11 

0.417391E+01 

0. 375000E+00 

0.551237E+01 

0. 136079E+01 

12 

0 . 360000Ef02 

0 . 607500Ef01 

0.585327E+02 

0.220362E+02 

13 

0.400000E+01 

0. 375000E+00 

0. 380878E+01 

0 . 234015E+00 

14 

0. 360000E+02 

0.607500Ef01 

0 . 327797Ef02 

0. 3794 16E-I-0 1 

15 

0 . 400000E+01 

0. 375000E+00 

0 . 533846E+01 

0. 136189E+01 

16 

0. 360000E+02 

0.607500E+01 

0.58S391Ef02 

0 . 220409E+02 

Note:  SOJOURN  refers  to  the  expected  sojourn  tiae  of  an  entity  to 
complete  the  entire  network.  WAIT(A)  refers  to  the  expected  waiting 
ttae  in  the  queue  at  the  fourth  node  in  the  network. 
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Appendix  D:  Tables  of  Results 


The  output  files  of  the  program  RESULTS  are  provided  as  tables  in 
this  appendix.  Tables  D.l  through  D.16  summarize  the  results  of  the  ten 
control  variates  against  the  response  variable  sojourn  time  at  the 
sixteen  design  points.  Similarly,  Tables  D.17  through  D.32  summarize 
the  results  of  the  ten  control  variates  against  the  response  variable  of 
the  probability  that  the  number  in  queue  exceeds  twice  the  expected 
number.  These  files  list  the  minimum,  maximum,  and  mean  values  of  the 
uncontrolled  response,  denoted  by  YBAR,  first.  Then  for  each  control 
variate  the  minimum,  maximum,  and  mean  values  of  the  controlled 
response,  denoted  by  YBAR(BHAT),  and  the  variance  ratio  are  listed. 

Also,  the  percentage  of  confidence  intervals  about  the  controlled 
response  that  cover  the  mean  value  of  YBAR  are  given  for  each  control 
variate. 
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Table  D.l.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  First  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERACE 


YBAR  0 . 369443E+01  0. 381331E+01  0.396763E+01 

USING  CONTROL  VARIABLE:  ROUTINC(l,3) 

YBAR(BHAT)  0.369317E+01  0.381367E*01  0.399705E+01  0.9300 

VARIANCE  RATIO  0.648974E+00  0.936955E+00  0.999987E+00 

USING  CONTROL  VARIABLE:  ROUTING* 1,4) 

YBAR(BHAT)  0.366235E+01  0.380847E+01  0.399742E+01  0.0850 

VARIANCE  RATIO  0.424270E+00  0 . 867965Ef00  0.999990E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.367848E+01  0.3811 38E+01  0.393465E+01  0.9200 

VARIANCE  RATIO  0.621534E+00  0.954747E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK* 2) 

YBAR(BHAT)  0.361416E+01  0.381202E+01  0.396768E+01  0.9250 

VARIANCE  RATIO  0.491398E+00  0.948515E+00  0 . 100000E+0 1 

USING  CONTROL  VARIABLE:  WORK* 3) 

YBAR(BHAT)  0.367779E+01  0.381282E+01  0.396515E+01  0.9350 

VARIANCE  RATIO  0.631606E+00  0.949113E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK* 4) 

YBAR(BHAT)  0.367833E+01  0.381140E+01  0.393472E+01  0.9200 

VARIANCE  RATIO  0.621459E+00  0.954740E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  SOJOURN(M/M/l ) 

YBAR(BHAT)  0.360603E+01  0.383068Ef01  0.412233E+01  0.5500 

VARIANCE  RATIO  0.125763E+00  0.418483E+00  0.817554E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.365523E+01  0.381422E+01  0.402782E+01  0.8300 

VARIANCE  RATIO  0.317536E+00  0.768369E+00  0 . 9999 78E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.361327E+01  0.382999E+01  0.412578E+01  0.5450 

VARIANCE  RATIO  0.859748E-01  0.384818E+00  0.774845E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.365208E+01  0.381421E+01  0.401822E+01  0.8100 

VARIANCE  RATIO  0.288632Ef00  0.751309Ef00  0.999933E+00 


Table  D.2.  Control  Variate  Results  Against  Sojourn  Tiae 
at  the  Second  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 253301E+02  0.292945E+02  0.340209E+02 

USING  CONTROL  VARIABLE:  R0UTING(1,3) 

YBAR(BHAT)  0.253703E+02  0.293262E+02  0.343614E+02  0.9350 

VARIANCE  RATIO  0 . 629330E+00  0.941911E+00  0.999992E+C0 

USINC  CONTROL  VARIABLE:  ROUTING* 1,4) 

YBAR(BHAT)  0.252341E+02  0.292181E+02  0.344292E+02  0.9550 

VARIANCE  RATIO  0.6U720E+00  0.906871E+00  0.999988E+00 

USING  CONTROL  VARIABLE:  WORK*l) 

YBAR(BHAT)  0 . 250874E+02  0 . 292620E+02  0.351617E+02  0.9400 

VARIANCE  RATIO  0.470587E+00  0.943858E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WORK* 2) 

YBAR(BHAT)  0.240332E+02  0.293179E+02  0.371701E+02  0.9350 

VARIANCE  RATIO  0.673789E4-00  0.946194Ef00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK* 3) 

YBAR*  BHAT)  0.253199E+02  0.293383Et02  0.339248E+02  0.9300 

VARIANCE  RATIO  0.608466E+00  0.950724E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  WORK* 4) 

YBAR* BHAT)  0.250668E+02  0.292583E+02  0.350784E+02  0.9400 

VARIANCE  RATIO  0.472088E+00  0.944131E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( N/M/ 1 ) 

YBAR*  BHAT)  0.260485E+02  0.307690E+02  0.381936E+02  0.8000 

VARIANCE  RATIO  0.152546E+00  0.697212E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR* BHAT)  0.247132E+02  0.305386E+02  0.365465E+02  0.6950 

VARIANCE  RATIO  0.157891E+00  0 . 644684E+00  0.999635E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR* BHAT)  0.260450E+02  0.307653E+02  0.383204E+02  0.8000 

VARIANCE  RATIO  0 . 150624E+00  0.698383E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  WAIT4(C/G/1> 

YBAR* BHAT)  0.247680E+02  0.305323E+02  0.365467E+02  0.7000 

VARIANCE  RATIO  0.159307E+00  0.641734E*00  0.995649E+00 


Table  D.3.  Control  Variate  Results  Against  Sojourn  Time 
at  the  Third  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERACE 

YBAR  0.477141E+01  0.521644E+01  0 .585939E+01 

USINC  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT)  0.473141E+01  0.521239E+01  0.580107E+01  0.94S0 

VARIANCE  RATIO  0.590194E+00  0.940776E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  R0UTINC(1,4) 

YBAR(BHAT)  0.478012E+01  0.521678E+01  0.590412E+01  0.9450 

VARIANCE  RATIO  0.537425E*00  0.925474E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  UORK(l) 

YBAR(BHAT)  0.483237E+01  0 . 521699E+01  0.582824E+01  0.9S00 

VARIANCE  RATIO  0.416803E+00  0 . 955932E+00  0.999999E+00 

USINC  CONTROL  VARIABLE:  UORK(2) 

YBAR(BHAT)  0.477196E+01  0.521219E+01  0.617023E+01  0.9500 

VARIANCE  RATIO  0.446203E+00  0.950714Et00  0.999980E+00 

USINC  CONTROL  VARIABLE:  UORK<3) 

YBAR(BHAT)  0.475023E+01  0.521933Et01  0.599340Ef0i  0.9450 

VARIANCE  RATIO  0.599295E+00  0.953930E+00  0.999999E+00 

USINC  CONTROL  VARIABLE:  WORK<4) 

YBAR(BHAT)  0.483239E+01  0.521705E+01  0.582838E+01  0.9500 

VARIANCE  RATIO  0.416409E+00  0.955908E+00  0.999997E+00 

USINC  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.463239E+01  0.528310E+01  0.655387E+01  0.5700 

VARIANCE  RATIO  0.125285E+00  0.468408E+00  0.890346E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/H/l) 

YBAR(BHAT)  0.455333E+01  0.527696E+01  0.645269E+01  0.5100 

VARIANCE  RATIO  0.100617E+00  0.386469E+00  0 . 791341E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( C/C/ 1 ) 

YBAR(BHAT)  0.451417E+01  0.527574E+01  0.645615Et01  0.3650 

VARIANCE  RATIO  0.620310E-01  0.250908E+00  0.667591E+00 

USINC  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR(BHAT)  0.439451E+01  0.526400E+01  0.635990E+01  0.4300 

VARIANCE  RATIO  0.699815E-01  0.316967E+00  0.703822E+00 


Table  D.4.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Fourth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERACE 

YEAR  0 . 371266E+02  0.457255E+02  0.5SS158E+02 

USING  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT)  0.378508E+02  0.457554E+02  0.548320E+02  0.9300 

VARIANCE  RATIO  0.486904E+00  0.942278E+00  0.999998E+00 

USINC  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.373728E+02  0.457098E+02  0.554014E+02  0.9350 

VARIANCE  RATIO  0.S31825Ef00  0.932978E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  U0RE(1) 

YBAR(BHAT)  0 . 375918E+02  0.458027E+02  0.623515E*02  0.9100 

VARIANCE  RATIO  0.652300E+00  0.947219E+00  0.999990E+00 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.344602E+02  0.456777E+02  0.573404E+02  0.9000 

VARIANCE  RATIO  0 .567247E+00  0 . 941918E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  UORK(3) 

YBAR(BHAT)  0.349138E+02  0.457327Et02  0.552863E+02  0.9150 

VARIANCE  RATIO  0.574240E+00  0.943468E+00  0.999986E+00 

USINC  CONTROL  VARIABLE:  W0RK(4) 

YBAR(BHAT)  0.376215E+02  0.457768E+02  0.618071E+02  0.9100 

VARIANCE  RATIO  0.651065E+00  0.947981E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.376422Ef02  0.498143E+02  0.691257Et02  0.6150 

VARIANCE  RATIO  0.160514E+00  0.662006E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.380585E+02  0.4921UE+02  0.671139E+02  0.6400 

VARIANCE  RATIO  0. 101392E*00  0.615236E+00  0.999963E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.376549E+02  0.498277E+02  0.694260E+02  0.5800 

VARIANCE  RATIO  0.731372E-01  0.600644E+00  0.999975E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/C/1) 

YBAR(BHAT)  0.382182E+02  0.492173E+02  0.668756E+02  0.6250 

VARIANCE  RATIO  0.971600E-01  0.610729E+00  0.999952E+00 


Table  D.5.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Fifth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERAGE 

YBAR  0 . 368183E+01  0.381449E+01  0.400045E+01 

USING  CONTROL  VARIABLE:  ROUTINGU.3) 

YBAR(BHAT)  0.363501E+01  0.381666E+01  0.401037E+01  0.9100 

VARIANCE  RATIO  0.407602E+00  0.869664E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  ROUTINC(l,4) 

YBAR(BHAT)  0.3SB449E+01  O.301251E+O1  0.406761E+01  0.9250 

VARIANCE  RATIO  0.388817E+00  0.88S916E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORR(l) 

YBAR(BHAT)  0 . 363304E+01  0.381513E+01  0.401350E+01  0.9350 

VARIANCE  RATIO  0.732037E+00  0.957379E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.367582E+01  0.381078E+01  0.402608E+01  0.9600 

VARIANCE  RATIO  0.769947E+00  0.957897E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORJU3) 

YBAR(BHAT)  0 . 358585E+01  0.381085E+01  0.398226E+01  0.9550 

VARIANCE  RATIO  0.707264E+00  0.954170E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.363303E+01  0.3B1514E+01  0.401370E+01  0.9350 

VARIANCE  RATIO  0 . 731882E+Q0  0.957361E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.361066E*01  0.384197E+01  0.414274E+01  0.4950 

VARIANCE  RATIO  0.139045E+00  0.377817E+00  0.881662E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.364710E+01  0.382226E+01  0.408382E+01  0.8150 

VARIANCE  RATIO  0 . 335014E+00  0.797228E+00  0.999920E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0. 361933E+-01  0.384037E+01  0.412813E+01  0.4750 

VARIANCE  RATIO  0.921910E-01  0 . 354609E+00  0.767348E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.364605E+01  0.382218E+01  0.407732E+01  0.8200 

VARIANCE  RATIO  0.347844E+00  0.783402E+00  0.999992E+00 


Table  D.6.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Sixth  Design  Point 


NAME  MINIMUM 

AVERAGE 

MAXIMUM 

COVERAGE 

‘  YBAR  0 . 248461E+02 

0.287763E+02 

0 . 332896E+02 

USING  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0.249434E+02  0.287970E+02 

VARIANCE  RATIO  0.6509B9E+00  0.941351E+00 

0 . 341 161E+02 
0.999999E+00 

0.9200 

'  USING  CONTROL  VARIABLE:  ROUTING! 1, 4 > 

’  YBAR(BHAT)  0.23S258E+02  0.287802E+02 

VARIANCE  RATIO  0.541487E+00  0.917009E+00 

0 . 348140E+02 
0.999993E+00 

0.9250 

USING  CONTROL  VARIABLE:  WORK(l) 
YBAR(BHAT)  0.2S1S09E+02 

VARIANCE  RATIO  0.663926E+00 

0 . 288086E+02 

0 . 959623E+00 

0 . 346363E+02 
0.100000E+01 

0.9450 

USING  CONTROL  VARIABLE:  WORK! 2) 
YBAR(BHAT)  0.245202E+02 

VARIANCE  RATIO  0.437341E+00 

0 . 287082E+02 

0 . 947478Ef00 

0 . 334179E+02 

0 . 100000E+0 1 

0.9400 

USING  CONTROL  VARIABLE:  WORK! 3) 
YBAR! BHAT)  0.243634E+02 

VARIANCE  RATIO  0.508923E+00 

* 

0 . 288033E+02 

0 .951356E+00 

0 . 370763E+02 

0 . 999982E+00 

0.9350 

USINC  CONTROL  VARIABLE:  WORK! 4) 
YBAR! BHAT)  0.251519E+02 

VARIANCE  RATIO  0.665197E*00 

0 . 288058E+02 
0.959802Ef00 

0.346073E+02 

0 . 100000E+01 

0.9450 

USINC  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR! BHAT)  0.257135E+02  0.302904E+02 

'  VARIANCE  RATIO  0.236069E+00  0.771186E+00 

0 . 382519E+02 
0.999962E+00 

0.8100 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR! BHAT)  0.241961E+02  0.300718E+02 

VARIANCE  RATIO  0.131538E*00  0.649533E+00 

0 . 380983E+02 

0 . 999407E+00 

0.7450 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBARIBHAT)  0.257132E+02  0.302892E+02 

VARIANCE  RATIO  0.245796E+00  0.771530E+00 

0.383264E+02 

0 . 999961E+00 

0.8100 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1> 

YBAR! BHAT)  0.243035E+02  0.300624E+02 

VARIANCE  RATIO  0.12S844E+00  0.644162E+00 

0 . 380028E+02 

0 . 999675E+00 

0.7200 

D.  7 


Table  D.7.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Seventh  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERACE 


YBAR  0.47311 8  E+0 1  0.521470E+01  0.582294E+01 

USING  CONTROL  VARIABLE:  ROUTING* 1,3) 

YBAR(BHAT)  0.475897E+01  0.521184E+01  0.592922E+01  0.9200 

VARIANCE  RATIO  0.572002E+00  0.927714E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  ROUTING(l,4) 

YBAR(BHAT)  0.469524E+01  0.521139E+01  0.585067E+01  0.8950 

VARIANCE  RATIO  0.608010E+00  0.914174E+00  0.999992E+-00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.469958E+01  0 . 521 1 12E+0 1  0.604592E+01  0.9300 

VARIANCE  RATIO  0.551904E+00  0.950072E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK* 2) 

YBAR(BHAT)  0.473291E+01  0.521587E+01  0.576332E+01  0.9250 

VARIANCE  RATIO  0.474959E+00  0.950610E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK* 3) 

YBAR(BHAT)  0.476462E+01  0.520558E+01  0.563329E+01  0.9050 

VARIANCE  RATIO  0.634175E+00  0.948817E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WORK* 4) 

YBAR(BHAT)  0.469965E+01  0.S21116E+01  0.604597E+01  0.9250 

VARIANCE  RATIO  0.551983Ef00  0.950057E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.464045E+01  0.527604E+01  0 . 620366E+01  0.6100 

VARIANCE  RATIO  0.968237E-01  0.507894E+00  0.951437E+00 

USING  CONTROL  VARIABLE:  WAIT4*M/N/1> 

YBAR(BHAT)  0.456468E+01  0.526524E+01  0.639802E+01  0.5450 

VARIANCE  RATIO  0.111213E+00  0.422013E*00  0.837698E+00 

USING  CONTROL  VARIABLE:  SOJOURN(G/G/l ) 

YBAR(BHAT)  0.450120E+01  0.526341E+01  0.640132E+01  0.4200 

VARIANCE  RATIO  0.863734E-01  0.279231E+00  0.658191E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.450152E+01  0.524690E+01  0.646256E+01  0.5100 

VARIANCE  RATIO  0.934624E-01  0.369244E+00  0.861650E+00 


D.  8 


Table  D.8.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Eight  Design  Point 

NAME  MINIMUM 

AVERACE 

MAXIMUM 

COVERAGE 

YBAR  0. 368068E+02 

0 . 451 1 13E+02 

0.556231E+02 

USING  CONTROL  VARIABLE:  ROUTING* 1,3) 

YBAR(BHAT)  0.367612E+02  0 . 45 1561E+02 

VARIANCE  RATIO  0.557175E+00  0.9399S7E+00 

0 . 568654E+02 

0 . 100000E+01 

0.9450 

USING  CONTROL  VARIABLE:  ROUTING* 1,4) 

YBAR(BHAT)  0.356995E+C2  0.451467E+02 

VARIANCE  RATIO  0.640790E+00  0.945581E+00 

0.553535E+02 

O.lOOOOOEfOl 

0.9550 

USING  CONTROL  VARIABLE:  WORK(l) 
YBAR(BHAT)  0.370118E+02 

VARIANCE  RATIO  0.615958E+00 

0.451506E+02 

0.951195E+00 

0.S63468E+02 

O.lOOOOOEfOl 

0.9450 

US  INC  CONTROL  VARIABLE:  WORK(2) 
YBAR(BHAT)  0.3691 3 3E+02 

VARIANCE  RATIO  0.492479E+00 

0.452114E+02 

0.947571E+00 

0.562391E+02 

0 . 999998E+00 

0.9300 

USING  CONTROL  VARIABLE:  WORK(3) 
Y8AR( BHAT)  0.369315Ef02 

i  VARIANCE  RATIO  0.552045E+00 

0 . 4S2633E+02 
0. 939136E+00 

0.572890E+02 

O.IOOOOOE+Ol 

0.9500 

US  INC  CONTROL  VARIABLE:  WORK(4) 
YBAR (BHAT)  0.370132E+02 

VARIANCE  RATIO  0.613340E+00 

0.451 325E+02 
0.952044E+00 

0.562732E+02 

0 . 100000E+01 

0.9450 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.363551E+02  0.495517E+02 

,  VARIANCE  RATIO  0.167355E+00  0.703501E+00 

0 . 659250E+02 
0. 998401E+00 

0.6550 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR (BHAT)  0.371424E+02  0.483798E+02 

VARIANCE  RATIO  0.105631E*00  0.620355E+00 

0 . 614S34E+02 
0 . 999962E+00 

0.6900 

USING  CONTROL  VARIABLE:  SOJOURN(C/C/l ) 

YBAR(BHAT)  0.363782E+02  0.494878E+02 

VARIANCE  RATIO  0 . 100811E+00  0.608829E+00 

0 .  635396Ef02 

0 . 999998E+00 

0.5900 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR (BHAT)  0.372278E+02  0.483759E+02 

VARIANCE  RATIO  0.104318Et00  0.6l7487Ef00 

0.623317E+02 

0 . 999963E+00 

0.6950 

.9 


Table  D.9.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Ninth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 30O592E+O1  0.398874E+01  0.415140E+01 

USING  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT)  0.382653E+01  0.398936E+01  0.415606E+01  0.9150 

VARIANCE  RATIO  0.650857E+00  0.949540E+00  0.999997E+00 

US  INC  CONTROL  VARIABLE:  ROUTINC(l,4) 

YBAR(BHAT)  0.375724E+01  0.399063E+01  0.420255E+01  0.7650 

VARIANCE  RATIO  0 . 37662SE+00  0.761644E+00  0 .998S50E+00 

USINC  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.380441E+01  0.398905E+01  0.417934E+01  0.8950 

VARIANCE  RATIO  0.513316E+00  0.938565E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORg(2) 

YBAR(BHAT)  0.379973E+01  0.398678E+01  0.416616E+01  0.9300 

VARIANCE  RATIO  0.587106E+00  0.948372E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK<3) 

YBAR(BHAT)  0.37S108E+01  0.398441E+01  0.418557E+01  0.9050 

VARIANCE  RATIO  0.634476E+00  0.943080E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.380442E+01  0.398906E+01  0.417932E+01  0.8950 

VARIANCE  RATIO  0.513602E+00  0.938575Ef00  0.999999E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.371078E+01  0.40l222EfQl  0.427863E+01  0.4650 

VARIANCE  RATIO  0.113236E+00  0.355449E+00  0.821358E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.370819E+01  0.399823E+01  0.424922E+01  0.6950 

VARIANCE  RATIO  0.205460E+00  0.662065E+00  0.999135E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.371430E+01  0.401247E+01  0.428437E+01  0.4400 

VARIANCE  RATIO  0.114329E+00  0.334087E+00  0.783583E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/C/1) 

YBAR(BHAT)  0.372064E+01  0.399890E+01  0.426272E+01  0.7100 

VARIANCE  RATIO  0.238847E+00  0 . 643083E+-00  0.997044E+00 


D.  10 


Table  D.10.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Tenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0.246467E+02 

USING  CONTROL  VARIABLE:  ROUTING( 
YBAR(BHAT)  0.247080E+02 

VARIANCE  RATIO  0.639244E+00 

USING  CONTROL  VARIABLE:  ROUTINC( 
YBAR(BHAT)  0.241610E+02 

VARIANCE  RATIO  0.293134Ef00 

USINC  CONTROL  VARIABLE:  UORK(l) 
YBAR(BHAT)  0.24S217E*02 

VARIANCE  RATIO  0.609402E+00 

USINC  CONTROL  VARIABLE:  WORK* 2) 
YBAR(BHAT)  0.236332E+02 

VARIANCE  RATIO  0.490976E+00 

USING  CONTROL  VARIABLE:  WORK(3) 
YBAR(BHAT)  0.247234Ef02 

VARIANCE  RATIO  0.635769E+00 

USING  CONTROL  VARIABLE:  WORK(4) 
YBAR(BHAT)  0.245203E+02 

VARIANCE  RATIO  0.609347E+00 


0 . 299427E+02 

0 . 353747E+02 

,3) 

0. 300389E+02 
0.951075E+00 

0 . 361594E+02 

0 . 999996E+00 

0.9100 

,*) 

0 . 299488E+02 

0 . 859462E+00 

0 . 357995E+02 

0 . 999998E+00 

0.8950 

0 . 299633E+02 
0.943036E+00 

0 . 363288E+02 
0.999994E+00 

0.8900 

0 . 299989E+02 

0 . 95 1623E+00 

0 . 357056E+02 

0 . 999996E+00 

0.9250 

0.299672E+02 

0 . 948150E+00 

0 . 363240E+02 

0 . 100000E+01 

0 . 8950 

0 . 299608E+02 

0 . 943254E+00 

0 . 363192E+02 

0 . 999999E+00 

0.8900 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.256705Ef02  0.314889E+02  0.411867E+02  0.7850 

VARIANCE  RATIO  0.273576E+00  0.708201E+00  0.100000E+01 


USING  CONTROL  VARIABLE:  WAIT4< M/M/1) 

YBAR(BHAT)  0.241567E+02  0.314347E+02  0.384832E+02  0.6850 

VARIANCE  RATIO  0 . 206467E+00  0.604245E+00  0.999674E+00 


USING  CONTROL  VARIABLE:  SOJOURN ( G/C/ 1 ) 

YBAR(BHAT)  0.256697E+02  0.314904E+02  0.413126E+02  0.7800 

VARIANCE  RATIO  0.277785E+00  0.709799E+00  0.999999E+00 


USING  CONTROL  VARIABLE:  WAIT4(C/G/1) 

YBAR(BHAT)  0.241370E+02  0.314233E+02  0.384162E+02  0.6600 

VARIANCE  RATIO  0.200296E+00  0.602308E+00  0.999712E+00 


D.  11 


Table  D.ll.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Eleventh  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERACE 


YBAR  0.490559E+01 

USING  CONTROL  VARIABLE:  ROUTING( 
YBAR(BHAT)  0.490509E+01 

VARIANCE  RATIO  0.572129E+00 

USING  CONTROL  VARIABLE:  ROUTINCC 
YBAR(BHAT)  0.491179E+01 

VARIANCE  RATIO  0.421668E+00 

USING  CONTROL  VARIABLE:  WORK(l) 
YBAR(BHAT)  0 . 490478E+01 

VARIANCE  RATIO  0.558922E+00 

USING  CONTROL  VARIABLE:  WORK(2) 
YBAR(BHAT)  0.4404S5E+01 

VARIANCE  RATIO  0.645532E*00 

USING  CONTROL  VARIABLE:  WORK(3) 
YBAR(BHAT)  0.490411E+01 

VARIANCE  RATIO  Q.516269E+00 

USING  CONTROL  VARIABLE:  WORK(4) 
YBAR(BHAT)  0.490473E+01 

VARIANCE  RATIO  0.557668E+00 


0 . 530441E+01 

0.566416E+01 

,3) 

0.530760E+01 

0 . 946071E+00 

0.566931E+01 

0.100000E+01 

0.9350 

,4) 

0.531109E+01 

0 . 873746E+00 

0.575136E+01 

0 . 100000E+01 

0.9050 

0.530375E+01 

0 . 946952E+00 

0.572119E+01 

0. 100000E+01 

0.9300 

0.529597E+01 

0 . 949085E+00 

0.582324E+01 

0.999978E+00 

0.9350 

0.530l30Ef01 

0.944746E+00 

0. 623082E+01 

0 . 999989E+00 

0.9350 

0.530378E+01 

0 . 946925E+00 

0.572085E*01 

0 . 1  OOOOOE-t-0 1 

0.9300 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.478901E+01  0 . 536703E+01  0.608927E+01  0.5800 

VARIANCE  RATIO  0.136079E+00  0.450433E+00  0.902546E+00 


USING  CONTROL  VARIABLE:  WAIT4< M/M/1) 

YBAR(BHAT)  0.466496E+01  0.535122E+01  0.598138E+01  0.5800 

VARIANCE  RATIO  0.936413E-01  0.415536E+Q0  0.979822E+00 


USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.461028E+01  0.535548E+01  0.608413E+01  0.4650 

VARIANCE  RATIO  0.653325E-01  0.266464E+00  0.702594E+00 


USING  CONTROL  VARIABLE:  WAIT4(C/C/1> 

YBAR(BHAT)  0.462090E+01  0.533501E+01  0.596814E+01  0.5250 

VARIANCE  RATIO  0.131066E+00  0.386874E*00  0.930679E+00 


D.  12 


Table  0.12.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Twelfth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 385956E+02  0.470301Ef02  0.590437E+02 

USING  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT)  0 . 38S054E*02  0.472231E+02  0.59573SE+02  0.9150 

VARIANCE  RATIO  0 . 678291E+00  0.949828E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTINCU.4) 

YBAR(BHAT)  0.375837E+02  0.468042E+02  0.589178E+02  0.9050 

VARIANCE  RATIO  0.54680SE+00  0.91 7021E+00  0  . 100000E+0 1 

USING  CONTROL  VARIABLE:  WORK( 1 ) 

YBAR(BHAT)  0.387353Ef02  0.471866E+02  0.688168E+02  0.9250 

VARIANCE  RATIO  0.600185E+00  0.949178E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  WORK<2) 

YBAR(BHAT)  0.385828E+02  0.470423E+02  0.596938E+02  0.9350 

VARIANCE  RATIO  0.615101E+00  0.953366E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.389498E*02  0.4711S9E+02  0.588525E+02  0.9200 

VARIANCE  RATIO  0.386973E+00  0.945840E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK<4) 

YBAR(BHAT)  0.387146E+02  0.471739E+02  0.679613E+02  0.9250 

VARIANCE  RATIO  0.595894Ef00  0.949470E+00  0.999991E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.393682E+02  0.505069E+02  0.698632E+02  0.6800 

VARIANCE  RATIO  0.135096E+00  0.689193E+00  0.999557E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.386547E+02  0.497819E+02  0.698244E+02  0.7250 

VARIANCE  RATIO  0.13S505E+00  0.671610E+00  0.998631E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( C/G/ 1 ) 

YBAR(BHAT)  0.393751E+02  0.502915E+02  0.699491E+02  0.6850 

VARIANCE  RATIO  0.129665E+00  0.670991E+00  0.999958E+00 

USINC  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.390337E+02  0.4977UE+02  0.699372E+02  0.7200 

VARIANCE  RATIO  0.143626E+00  0.672896E+00  0 . 998564E+00 


0.13 


Table  D.13.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Thirteenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 362597E+01  0.381095E+01  0.398981E+01 

USING  CONTROL  VARIABLE:  R0UTING(1,3) 

YBAR(BHAT)  0.362032E+01  0 . 38091 1E+01  0.400185E+01  0.8800 

VARIANCE  RATIO  0.335364E+00  0.877644E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.359593E+01  0.381070E+01  0.406907E+01  0.8000 

VARIANCE  RATIO  0.212SS8E+00  0.776S82E+00  0 . 100000E+01 

USING  CONTROL  VARIABLE:  WORR(l) 

YBAR(BHAT)  0.362653E*01  0.3810S4Et01  0.398310E+01  0.8950 

VARIANCE  RATIO  0.487505E+00  0.94S541E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORR(2) 

YBAR(BHAT)  0.362480E+01  0.381052E+01  0.398868E+01  0.9100 

VARIANCE  RATIO  0.385303E+00  0 . 9491 18E+00  O.lOOOOOEtOl 

USING  CONTROL  VARIABLE:  WORK<3) 

YBAR(BHAT)  0.361801E+01  0.380606E+01  0.399455E+01  0.9150 

VARIANCE  RATIO  0.476760E+00  0.942032E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORSU4) 

YBAR(BHAT)  0.362652E+01  0.3810S4E+01  0.398302E+01  0.8950 

VARIANCE  RATIO  0.487568E+00  0.945545E+00  O.lOOOOOEtOl 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.354936E+01  0.383049E+01  0.414142E+01  0.4300 

VARIANCE  RATIO  0.611835E-01  0.331030E+00  0.881636E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.358877E+01  0.381797E+01  0.403301E+01  0.7650 

VARIANCE  RATIO  0.179257E+00  0.712198E+00  0.999670E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.354646E+01  0.382954E+01  0 . 41 361 2E+01  0.4350 

VARIANCE  RATIO  0.643630E-01  0.315746E+00  0.865978E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/C/1) 

YBAR(BHAT)  0.356530E+01  0.381773E+01  0.404205E+01  0.7400 

VARIANCE  RATIO  0 . 19i999E*00  0.699534E+00  0.999399E+00 


D.  14 


Table  D.14.  Control  Variate  Results  Against  Sojourn  Time 
at  the  Fourteenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YEAR  0.250717E+02  0.292389E+02  0.348594E+02 

USING  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0.244769E+02  0.292682E*02  0.348S73E+02  0.9450 

VARIANCE  RATIO  0.417618E+00  0.928636E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR(BHAT)  0.245819E+02  0.293296E+02  0.367279E+02  0.9000 

VARIANCE  RATIO  0 .443243E4-00  0.858234Et00  0.999979E+00 

USING  CONTROL  VARIABLE:  WORR(l) 

YBAR(BHAT)  0.252217Ef02  0.292557E+02  Q.361015E+02  0.9450 

VARIANCE  RATIO  0.506835E+00  0.952300E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR(BHAT)  0.231941E+02  0.292391E+02  0.376504E+02  0.9300 

VARIANCE  RATIO  0.702847E+00  0.945412E+00  0 . 999991E+00 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR(BHAT)  0.257546E+02  0.293223E+02  0.348498E+02  0.9350 

VARIANCE  RATIO  0.562071E+00  0.948348E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WORK! 4) 

YBAR! BHAT '  0.252241E+02  0.292533E+02  0.360774E+02  0.9450 

VARIANCE  RATIO  0.510453E+00  0.952444E+00  0 . 999997E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR!  BHAT)  0 . 246920E4-02  0.3075 1 7E+02  0.384403E+02  0.8150 

VARIANCE  RATIO  0.173840E+00  0.748431E+00  0.999981E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR! BHAT)  0.224575E+02  0.306417E+02  0.385982E+02  0.7450 

VARIANCE  RATIO  0.205944E+00  0.632572E+00  0.995176E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR! BHAT)  0 . 244717E+02  0.307440E+02  0.383845E+02  0.8150 

VARIANCE  RATIO  0.156665E+00  0.750094E+00  0.999974E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR! BHAT)  0.219894E+02  0.306433E+02  0.385125E+02  0.7450 

VARIANCE  RATIO  0 . 195685E+00  0.626817E+00  0.993579E+00 


D.  15 


Table  0.15.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Fifteenth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERAGE 


YBAR  0.470558E+Q1 

USING  CONTROL  VARIABLE:  R0UTING( 
YBAR(BHAT)  0.470723E+01 

VARIANCE  RATIO  0.539615E+00 

USINC  CONTROL  VARIABLE:  ROUTING( 
YBAR(BHAT)  0.464156E+01 

VARIANCE  RATIO  0 .435000E4-00 

USING  CONTROL  VARIABLE:  WORK(l) 
YBAR(BHAT)  0.460887E+01 

VARIANCE  RATIO  0 .491643E+00 

USINC  CONTROL  VARIABLE:  W0RK(2) 
YBAR(BHAT)  0.445197E+01 

VARIANCE  RATIO  0.69S184E+00 

USING  CONTROL  VARIABLE:  WORK<3) 
YBAR(BHAT)  0.463500E+01 

VARIANCE  RATIO  0.615519E+00 

USINC  CONTROL  VARIABLE:  W0RK(4) 
YBAR(BHAT)  0.460917E+01 

VARIANCE  RATIO  0.491485E+00 


0.512254E+01 

0.S453S2E+01 

,3) 

0.51 25  0  7  E+ 0 1 
0.935454E+00 

0.545784E+01 

0 . 100000E+01 

0.9300 

,4) 

0.512207E+01 

0 . 899192E+00 

0 .551342E+01 

0 . 999996E+00 

0.9250 

0.511 759E+0 1 

0 . 939493E+00 

0.548528E+01 

0 . 100000E+01 

0.9150 

0.511944E+01 

0 ,957703Ef00 

0 . 548606E+01 

0 . 999999E+00 

0.9250 

0 . 5 12162E+0 1 
0.947813E+00 

0.546696E+01 

0 . 100000E+01 

0.9250 

0.511760E+01 

0 . 939497E+00 

0.548530E+01 

0. 100000E+01 

0.9150 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.446228E+01  0.517711E+01  0.578296E+01  0.5650 

VARIANCE  RATIO  0.123806E+00  0.448358E+00  0.937828E+00 


USING  CONTROL  VARIABLE:  WAIT4(M/M/1> 

YBAR(BHAT)  0.455079E+01  0.516825E+01  0.569372E+01  0.5150 

VARIANCE  RATIO  0 . 123214E*00  0.404259E+00  0.765727E+00 


USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.442982E+01  0.517591E+01  0.579209E+01  0.3850 

VARIANCE  RATIO  0.328271E-01  0.253222E+00  0.724730E+00 


USING  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR(BHAT)  0.450089E+01  0.516114E+01  0.572054E+01  0.4950 

VARIANCE  RATIO  0.938772E-01  0.379051E+00  0.778870E+00 


Table  D.16.  Control  Variate  Results  Against  Sojourn  Tine 
at  the  Sixteenth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERAGE 


YBAR  0 . 368060E+02  0.456687E+02  0.572155E+02 

USING  CONTROL  VARIABLE:  R0UTINC(1,3) 

YBAR(BHAT)  0.367545E+02  0.455973E*02  0.581674E+02  0.9150 

VARIANCE  RATIO  0.569230E+00  0.941763E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTING(l,4) 

YBAR(BHAT)  0.361866E+02  0.455062E+02  0.576799E+02  0.9050 

VARIANCE  RATIO  0.459856E+00  0.928508E+00  0.999994E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.371077E*Q2  0.457249E+02  0.600144E+02  0.89S0 

VARIANCE  RATIO  0 . 608909E+00  0.944781E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.370771E+02  0.457057E+02  0.576533E+02  0.9150 

VARIANCE  RATIO  0.646092E+00  0.944484Ef00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORR(3) 

YBAR(BHAT)  0.370829E+02  0.456703E*02  0.576092E+02  0.9050 

VARIANCE  RATIO  0.525826E+00  0.951012E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.371067E+02  0.457137E+02  0.598441E+02  0.8900 

VARIANCE  RATIO  0.608657E+00  0.945294E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.374931E+02  0.491399E+02  0.681191E+02  0.6550 

VARIANCE  RATIO  0.214510E+00  0.725930E+00  0.999967E+00 

USINC  CONTROL  VARIABLE:  WAIT4< M/M/1) 

YBAR(BHAT)  0.322618E+02  0.484143E+02  0.565718E+02  0.6500 

VARIANCE  RATIO  0.200033E-03  0.623759E+00  0.999869E+00 

USING  CONTROL  VARIABLE:  SOJOURN(C/C/l ) 

YBAR(BHAT)  0.374865E+02  0.491428E+02  0 . 688984E+02  0.6300 

VARIANCE  RATIO  0.700837E-02  0.655286E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/G/1) 

YBAR(BHAT)  0.318773E+02  0.483783E+02  0.663283E+02  0.6650 

VARIANCE  RATIO  0.200033E-03  0.624382E+00  0.999870E+00 


Table  D.17.  Control  Variate  Results  Against  Fourth  Hode  Quantile 

at  the  First  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERACE 


Y8AR  0. 108982E+00  0.121346E+00  0. 133839E*00 

USING  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0 . 107888E+00  0.1Z1325E+00  0 . 1 33859E+00  0.90S0 

VARIANCE  RATIO  0.S84686E+00  0.939366E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR(BHAT)  0.10694SE+00  0.121218E+00  0.133643E+00  0.8900 

VARIANCE  RATIO  0.332811E+00  0.909225E*00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.108926E+00  0.121482E+00  0.134308E+00  0.9000 

VARIANCE  RATIO  0.697780E+00  0.9S0873E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR(BHAT)  0.108073E+00  0.121293E+00  0.133982E+00  0.9100 

VARIANCE  RATIO  0 . 621506E+00  0.946486E+00  0.999983E+00 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR(BHAT)  0 . 10627SE+00  0.121373E+00  0.134094E+00  0.8950 

VARIANCE  RATIO  0.565127E+00  0.951841E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 4) 

YBAR(BHAT)  0.108926E+00  0.121483E+00  0.134307E+00  0.9000 

VARIANCE  RATIO  0.698097E+00  0.950868E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.106004Ef00  0.122063Ef00  0.137733Ef00  0.7500 

VARIANCE  RAixO  0.282425E+00  0.704738E+00  0.999972E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR! BHAT)  0.102571E+00  0.121814E+00  0.141884E+00  0.5700 

VARIANCE  RATIO  0.1S9547Ef00  0.443010E+00  0.837345E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR! BHAT)  0.106520E+00  0.121954E+00  0.138474E+00  0.7700 

VARIANCE  RATIO  0.283256E+00  0.720013E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WAIT4!C/G/1) 

YBAR! BHAT)  0.102915E+00  0.121783E+00  0.143014E+00  0.5600 

VARIANCE  RATIO  0.140546E+00  0.421196Ef00  0.832287E+00 


Table  D.18.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Second  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 

YBAR  0.800609E-01  0.122987E+00  0.183393E+00 

USING  CONTROL  VARIABLE:  R0UTING(1,3) 

YBAR(BHAT)  0.796589E-01  0.122874E+00  0.188708E+00  0.9400 

VARIANCE  RATIO  0.665956E+00  0.944736E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.690142E-01  0.122216E+00  0.187343E+00  0.9300 

VARIANCE  RATIO  0 .505389E+00  0.918495E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORR(l) 

YBAR(BHAT)  0.713SS5E-01  0.123318E+00  0.183295E+00  0.9300 

VARIANCE  RATIO  0 .548735E+00  0.946070E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORR<2) 

YBAR(BHAT)  0.708987E-01  0.122474E+00  0.188564E+00  0.9200 

VARIANCE  RATIO  0.520072E+00  0.944244E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.810515E-01  0.123990E+00  0.187336E+00  0.9350 

VARIANCE  RATIO  0.711601E+00  0.947854E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORK<4) 

YBAR(BHAT)  0.715930E-01  0.123326E+00  0.183302E+00  0.9300 

VARIANCE  RATIO  0.550165E+00  0.945971E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.767694E-01  0.127122E+00  0 . 197354E+00  0.9150 

VARIANCE  RATIO  0.381408E+00  0.918162E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.803986E-01  0 . 1 341 37E+00  0.213651E+00  0.7650 

VARIANCE  RATIO  0.247330E+00  0.723061E+00  0.999982E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.766782E-01  0.126980E+00  0.196664E+00  0.9150 

VARIANCE  RATIO  0 . 381622E+00  0.919540E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR(BHAT)  0.800893E-01  0.134095E+00  0.215977E+00  0.7600 

VARIANCE  RATIO  0 . 25 1796E+00  0.721281E+00  0.999996E+00 


Table  D.19.  Control  Variate  Results  Against  Fourth  Node  Quantile 

at  the  Third  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 1 17299E+00  0.141142E+00  0.165164E+00 

USING  CONTROL  VARIABLE:  ROUT INC (1,3) 

YBAR(BHAT)  0.114428E+00  0 . 140852E+00  0.167499E+00  0.9400 

VARIANCE  RATIO  0.647251E+00  0.946466E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  ROUTING(l,4) 

YBAR(BHAT)  0.117293E+00  0.14U51E+00  0.168151E+00  0.9300 

VARIANCE  RATIO  0.641702E+00  0.942667Ef00  0.10000QE+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.119165E+00  0.141217E+00  0.167973E+00  0.95S0 

VARIANCE  RATIO  0.464607E+00  0.9S170SE+00  0.999995E+00 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.116956E+00  0.140959E+00  0.169834E+00  0.9400 

VARIANCE  RATIO  0.532538E+00  0.951487E+00  0.999959E+00 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.115201E+00  0.141266E+00  0.167604E+00  0.9300 

VARIANCE  RATIO  0.555494E+00  0.950901E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  W0RK<4) 

YBAR(BHAT)  0.119136E+00  0.141222E+00  0.167956E+00  0.9550 

VARIANCE  RATIO  0.464179E+00  0.951710E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0 . 109585E+00  0.144851E*00  0 . 193024E+00  0.6600 

VARIANCE  RATIO  0.152195E+00  0.579707E+00  0.992372E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.939711E-01  0.145167E+00  0.193472E+00  0.4350 

VARIANCE  RATIO  0.809851E-01  0.268138E+00  0.677336E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.933476E-01  0.144652E+00  0.191033E+00  0.3650 

VARIANCE  RATIO  0.796989E-01  0.262950E+00  0.661056E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/C/i) 

YBAR(BHAT)  0.835182E-01  0.144021E+00  0.187470E+00  0.3750 

VARIANCE  RATIO  0.822466E-01  0 . 232744E+00  0.625153E+00 


Table  D.20.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Fourth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 738233E-02  0 . 613441E-01  0.141887E+00 

USING  CONTROL  VARIABLE:  ROUT INC (1,3) 

YBAR(BHAT)  0.758367E-02  0.614113E-01  0.136472E+00  0.8S00 

VARIANCE  RATIO  0.637302E+00  0.951017E+00  0.999995E+00 

USING  CONTROL  VARIABLE:  ROUTINC(l,4> 

YBAR(BHAT)  0.732101E-02  0.607389E-01  0.147175E+00  0.8650 

VARIANCE  RATIO  0.623753E+00  0.956300E+00  0.99998SE+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.486476E-02  0.615487E-01  0 . 163884E+00  0.8550 

VARIANCE  RATIO  0 . 339997E+00  0.944999E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORR(2) 

YBAR(BHAT)  -0.128170E-01  0.613136B-01  0.144628E+00  0.8400 

VARIANCE  RATIO  0.485296EfOO  0.946526E+00  0.999993E+00 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.530243E-02  0.617015E-01  0.148l24Et00  0.8450 

VARIANCE  RATIO  0.617724E+00  0.940743E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.486544E-02  0.614900E-01  0.163384E+00  0.8550 

VARIANCE  RATIO  0 . 339796E+00  0.945119E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.826074E-02  0.754763E-01  0.201013E+00  0.7350 

VARIANCE  RATIO  0.421055E-01  0.837181E+00  0.999949E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.747168E-02  0.771357E-01  0.201252E+00  0.6850 

VARIANCE  RATIO  0.162178E-01  0.751506E+00  0 . 999986E+00 

USING  CONTROL  VARIABLE:  SOJOURN(C/C/l ) 

YBAR(BHAT)  0.829340E-02  0.777838E-01  0.200373E+00  0.6850 

VARIANCE  RATIO  0.356739E-01  0.774744E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.769022E-02  0.771427E-01  0.197136E+00  0.6750 

VARIANCE  RATIO  0.162676E-01  0.747599E+00  0.100000E+01 


Table  D.21.  Control  Variate  Results  Against  Fourth  Node  Quantile 

at  the  Fifth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERACE 


YBAR  0. 111325E+00  0.121418E+00  0.132414E+00 

USING  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0.111788E+00  0.121434E+00  C.131716E+00  0.9100 

VARIANCE  RATIO  0.S81385E+00  0.947085E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR(BHAT)  0.105001E+00  0.121351E+00  0.132163E+00  0.8850 

VARIANCE  RATIO  0.S78932E+00  0.903681E+00  0.999983E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR( BHAT)  0 . 109S  UEfOO  0.121498E+00  0.133965E+00  0.9100 

VARIANCE  RATIO  0.702469E+00  0.949806E+00  0.999995E+00 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR (BHAT)  0.107856E+00  0.121264E+00  0.134142E+00  0.9050 

VARIANCE  RATIO  0.673564E4-00  0.953347E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR! BHAT)  0.105744E+00  0.121409E+00  0.134251E+00  0.9050 

VARIANCE  RATIO  0.627567E+00  0.95U25E+00  Q.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 4) 

YBAR! BHAT)  0.109515E+00  0.121498E+00  0.133964E+00  0.9100 

VARIANCE  RATIO  0.702539E+00  0.949797E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR! BHAT)  0.110182E+00  0.122408E+00  0.136964E+00  0.7800 

VARIANCE  RATIO  0.303391E+00  0.779417E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR (BHAT)  0.105685E+00  0.122206E+00  0.140539E+00  0.5450 

VARIANCE  RATIO  0.130080E+00  0.445908E+00  0.925947E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( C/C/ 1 ) 

YBAR! BHAT)  0.110200E+00  0.122268E+00  0.137154E+00  0.7950 

VARIANCE  RATIO  0.323205E+00  0.796578E+00  0.999815E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR! BHAT)  0.105442E+00  0.122140E+00  0.139733E+00  0.5100 

VARIANCE  RATIO  0.127179E+00  0.418037E+00  0.908348E+00 


D.  22 


Table  D.22.  Control  Variate  Results  Against  Fourth  Node  Quantile 

at  the  Sixth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0.622885E-01  0.106684E+00  0.168433E+00 

USING  CONTROL  VARIABLE:  R0UTING( 1,3) 

YBAR(BHAT)  0.620325E-01  0.108875EfQ0  0.17Q645E+00  0.92S0 

VARIANCE  RATIO  0.513468E+00  0.9506S7E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  R0UTINC(1,4> 

YBAR(BHAT)  0.676376E-01  0.108407E+00  0.169864E+00  0.9350 

VARIANCE  RATIO  0.583689E+00  0.939187E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.594641E-01  0.109773E+00  0 . 172394E4-00  0.9300 

VARIANCE  RATIO  0.612138E+00  0.953165E+00  0.999997E+00 

US  INC  CONTROL  VARIABLE:  WORK<2> 

YBAR(BHAT)  0.622200E-01  0.108051E+00  0.172440E+00  0.9250 

VARIANCE  RATIO  0.593514E+00  0.949890E+00  0.999988E+00 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.S87167E-01  0.108695E+00  0.180395E+00  0.9300 

VARIANCE  RATIO  0.468669E+00  0.948411E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.595041E-01  0.109778E+00  0.172424E+00  0.9300 

VARIANCE  RATIO  0.612476E+00  0.953100E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.580289E-01  0.110269Ef00  0.177221E+00  0.9150 

VARIANCE  RATIO  0.454583E+00  0.940221E+00  0.999979E+00 

USING  CONTROL  VARIABLE:  WAIT4( M/M/1) 

YBAR(BHAT)  0.651024E-01  0.120729E+00  0.223535E+00  0.8000 

VARIANCE  RATIO  Q.264705E+00  0.735405E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( G/C/ 1 ) 

YBAR(BHAT)  0.579146E-01  0.110227E+00  0.176488E+00  0.9150 

VARIANCE  RATIO  0.470616E+00  0.940779E+00  0.999979E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/C/i) 

YBAR(BHAT)  0.647745E-01  0.120575E+00  0.222433E+00  0.8000 

VARIANCE  RATIO  0.264302E+00  0.731573E+00  0.100000E+01 


D.23 


Table  D.23.  Control  Variate  Results  Against  Fourth  Hode  Quantile 
at  the  Seventh  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0. 1 1 2569E+00  0.140501E+00  0 . 185 1S6E+00 

USING  CONTROL  VARIABLE:  R0UTINC<1,3) 

YBAR(BHAT)  0. 113335E+00  0.140270Ef00  0.190232E+00  0.9400 

VARIANCE  RATIO  0.547357E+00  0.942611E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.109996E+00  0.140422E+00  0.186807E+00  0.9300 

VARIANCE  RATIO  0.S70730E+00  0.934657E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORR(l) 

YBAR(BHAT)  0 . 1 12687E+00  0.140S80E+00  0.201817E+00  0.9250 

VARIANCE  RATIO  0.607985E+00  0.944603E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.112672E+00  0.140510E+00  0.183614E+00  0.9300 

VARIANCE  RATIO  0.6U688E+00  0.948069E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.110439E+00  0.140202E+00  0.171872E+00  0.9150 

VARIANCE  RATIO  0.660412E+00  0.951762E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  W0RK(4) 

YBAR(BHAT)  0.112691E+00  0.140583E+00  0.201832E+00  0.9250 

VARIANCE  RATIO  0.608191E+00  0.944580E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  SOJOURN(M/M/1 ) 

YBAR(BHAT)  0. 108940E+-00  0.143605E+00  0.204380E+00  0.7500 

VARIANCE  RATIO  0.234856E+00  0.669272E+00  0.997503E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/l) 

YBAR(BHAT)  0.987313E-01  0.143796E+00  0.217044E+00  0.4850 

VARIANCE  RATIO  0.848460E-01  0.270045Ef00  0.888704E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.984601E-01  0.143404E+00  0.216561E+00  0.4850 

VARIANCE  RATIO  0.797107E-01  0.309553E+00  0.768329E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/G/1) 

YBAR(BHAT)  0.943801E-01  0.142356E+00  0.221662E+00  0.4350 

VARIANCE  RATIO  0.624765E-01  0.243331E+00  0.717937E+00 


D.24 


Table  0.24.  Control  Variate  Results  Against  Fourth  Node  Quantile 

at  the  Eight  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERAGE 


YBAR  0 . 145248E-01  0.645567E-01  0.177051E+00 

USING  CONTROL  VARIABLE:  R0UTING(1,3) 

YBAR(BHAT)  0.152320E-01  0.641725E-01  0.177135E+00  0.8700 

VARIANCE  RATIO  0.709700E+00  0.943849E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR(BHAT)  0.137351E-01  0.651024E-01  0.180345E+00  0.8650 

VARIANCE  RATIO  0.570944E+00  0.933924E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  -0.429269E-02  0.647223E-01  0.178368E+00  0.8650 

VARIANCE  RATIO  0.702579E+00  0.956512E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR(BHAT)  0.955922E-02  0.649550E-01  0.183485E+00  0.8600 

VARIANCE  RATIO  0 . 38815 1E+00  0.950948E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR! BHAT)  0.145098E-01  0.6S6369E-01  0 . 170770E+00  0.8650 

VARIANCE  RATIO  0.389760E+00  0.948317E*00  0.999971E+00 

USING  CONTROL  VARIABLE:  WORK! 4) 

YBAR! BHAT)  -0 . 441317E-02  0.646924E-01  0.178622E+00  0.B65C 

VARIANCE  RATIO  0.711Q96E+00  0.9S6903E+00  0.999995E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR! BHAT)  0.123770E-01  0.800651E-01  0.190834E+00  0.7550 

VARIANCE  RATIO  0.256196E+00  0.863874E+00  0.999977E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR! BHAT)  0.139049E-01  0.804204E-01  0.191781E+00  0.7200 

VARIANCE  RATIO  0.909681E-01  0.747265E+00  0.999910E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/C/ 1 ) 

YBAR (BHAT)  0.135561E-01  0.829255E-01  0.198735E+00  0.6700 

VARIANCE  RATIO  0.196644E+00  0.775468E+00  0.999973E+00 

USING  CONTROL  VARIABLE:  WAIT4( C/C/1) 

YBAR! BHAT)  0.139048E-01  0.804041E-01  0.192440E+00  0.7150 

VARIANCE  RATIO  0.116133E+00  0.746031E+00  0.100000E+01 
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Table  D.2S.  Control  Variate  Results  Against  Fourth  Mode  Quantile 

at  the  Ninth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 

YBAR  0.  Ill 008E+-00  0.122395E*00  0.132180E+00 

USING  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT)  0.109439E+00  0 . 122366E+-00  0.1330S2E+00  0.9300 

VARIANCE  RATIO  0.497627E+00  0.942552E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTINC(l,4) 

YBAR(BHAT)  0.109022E+00  0.122632E+00  0.134885E+00  0.8500 

VARIANCE  RATIO  0.433837E+00  0.838388E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.110767E+00  0.122497Ef00  0.134761Ef00  0.9100 

VARIANCE  RATIO  0.715594E+00  0.948371E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORJU2) 

YBAR(BHAT)  0.109774E+00  0.122348E+00  0.131130E+00  0.9400 

VARIANCE  RATIO  0.471377E+00  0 . 95 1846E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK(3) 

YBAR(BHAT)  0.U1225E+00  0.122125E+00  0.132695E+00  0.9250 

VARIANCE  RATIO  0.687864E+00  0.941987E+00  0.999993E+00 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.110766E+00  0.122497E+00  0.134753E+00  0.9100 

VARIANCE  RATIO  0. 715583E+-00  0.948372E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.106708E+00  0.123S62E*G0  0 . i 37 i**on+00  0.6900 

VARIANCE  RATIO  0.227410E+00  0.611888E+00  0.999926E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.101738E+00  0 . 1231 72E+00  0.139209E+00  0.4950 

VARIANCE  RATIO  0.104222Et00  0.395685E+00  0.840172E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.106979E+00  0.123520E+00  0.135986E+00  0.6900 

VARIANCE  RATIO  0.223015E+00  0.624305E+00  0.999225E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR(BHAT)  0.101454E+00  0.123210E+00  0.138830E+00  0.4650 

VARIANCE  RATIO  0.796724E-01  0 . 374487E+00  0.844993E+00 
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Table  D.26.  Control  Variate  Results  Against  Fourth  Node  Quantile 

at  the  Tenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 685479E-01  0.112705Ef00  0.177696E+00 

USING  CONTROL  VARIABLE:  ROUTINC(l,3> 

YBAR(BHAT)  0.645306E-01  0 . 1 1257SE+00  0 . 1 79238E4-00  0.9350 

VARIANCE  RATIO  0.639193E+00  0.953142E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR(BHAT)  0.657353E-01  0.112503E+00  0.175888E+00  0.9200 

VARIANCE  RATIO  0 .473448E+00  0.913959E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.648471E-01  0.114017E+00  0.177704E+00  0.9300 

VARIANCE  RATIO  0.669707E+00  0.951018E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR(BHAT)  0.665513E-01  0.112739E+00  0.177008E+00  0.9400 

VARIANCE  RATIO  0.703021E+00  0.955436E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR(BHAT)  0.546768E-01  0.112561E+00  0.182294E+00  0.9250 

VARIANCE  RATIO  0.678082E+00  0.953411E+00  0.100000E+01 

US  INC  CONTROL  VARIABLE:  WORK! 4) 

YBAR(BHAT)  0.647503E-01  0.114022E+00  0.177704E+00  0.9300 

VARIANCE  RATIO  0.668615E+00  0.951008E+00  0.999991Et00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.668401E-01  0 . 1 1695  7E+00  0.206801E+00  0.9300 

VARIANCE  RATIO  0.358874E+00  0.899939E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/i) 

YBAR(BHAT)  0.561665E-01  0.125670E+00  0.225853E+00  0.7600 

VARIANCE  RATIO  0.161109E+00  0 . 713087E+00  0.999976E+00 

USING  CONTROL  VARIABLE:  SOJOURN! G/G/l ) 

YBAR(BHAT)  0.668441E-01  0.116895E+00  0.206502E+00  0.9300 

VARIANCE  RATIO  0.364805E+00  0.901445E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WAIT4(G/G/1) 

YBAR(BHAT)  0.535490E-01  0.125568E+00  0.224347E+00  0.7500 

VARIANCE  RATIO  0.160630E+00  0.711324E+00  0.999846E+00 


Table  D.27.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Eleventh  Design  Point 


NAME 

MINIMUM 

AVERAGE 

MAXIMUM 

COVERACE 

YBAR 

0 . 109092E+00 

0.133791E+Q0 

0 . 1S63S6E+00 

USING  CONTROL  VARIABLE:  ROUTING(l,3) 

YBAR(BHAT) 

0 . 109037E+00 

0 . 133896E+00 

0 . 153300E+0G 

0.9650 

VARIANCE  RATIO 

0 . 717994E+00 

0 . 946812E+00 

0 . 999999E+00 

USING  CONTROL  VARIABLE:  ROUT INC (1,4) 

YBAR(BHAT) 

0 . 109303E+00 

0 . 134058E+00 

0 . 1S6681E+00 

0.9600 

VARIANCE  RATIO 

0 .583121E+00 

0 . 933723E+00 

0.999992E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT) 

0 . 105861E+00 

0 . 133454E+00 

0 . 160838E+00 

0.9350 

VARIANCE  RATIO 

0.595487E+00 

0 . 944692E+00 

0 . 999999E+00 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT) 

0 . 858128E-01 

0 . 1 33426E+00 

0 . 16415  3E+00 

0.9550 

VARIANCE  RATIO 

0.565806Ef00 

0 . 946036E+00 

0 . 999998E+00 

USING  CONTROL  VARIABLE:  WORR(3) 

YBAR(BHAT) 

0 . 1101 8  3E+-00 

0 . 133773E+00 

0 . 194814E+00 

0.9300 

VARIANCE  RATIO 

0.583655E+00 

0. 94619DE+00 

0. 100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT) 

0 . 105856E+00 

0. 133455E+00 

0 . 160815E+00 

0.9350 

VARIANCE  RATIO 

0 . 595 146Ef00 

0 . 944685E+00 

0 . 100000E  +  01 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT) 

0 . 103230E+-00 

0. 137180E+00 

0. 179016E+00 

0.7100 

VARIANCE  RATIO 

0 . 15  3193E+00 

Q.610513E+00 

0 . 995094E+00 

US  INC  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT) 

0 . 918101 E— 0 1 

0 . 1 37244E+00 

0 . 180314E+00 

O 

-J 

o 

o 

VARIANCE  RATIO 

0 . 74085  IE-01 

0 . 271448E+00 

0 . 791968E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT) 

0.906200E-01 

0 . 137089E+00 

0 . 1 79633E+00 

0.4850 

VARIANCE  RATIO 

0 . 659187E-01 

0.271682E+00 

0 . 65061 1E+0G 

USING  CONTROL  VARIABLE:  WAIT4(G/C/1) 

YBAR(BHAT) 

0 . 878249E-01 

0 . 1 35907E+00 

0. 178790E+00 

0.4300 

VARIANCE  RATIO 

0.654438E-01 

0.2387S1E+00 

0 . 696407E+00 
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Table  D.28.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Twelfth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 

YBAR  0 . 18301 3E-01  0.654637E-01  0.151707E+00 

USINC  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0.186208E-01  0.660864E-01  0.162674E+00  0.8750 

VARIANCE  RATIO  0.438063E+00  0.946059E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  ROUT INC (1,4) 

YBAR(BHAT)  0.180504E-01  0.647238E-01  0.154406E+00  0.8600 

VARIANCE  RATIO  0.631840E+00  0.934873Ef00  0999966E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.136727E-01  0.662149E-01  0.2081S2E+00  0.8800 

VARIANCE  RATIO  0.479767E+00  0.947534E+00  0.999994E+00 

USING  CONTROL  VARIABLE:  W0RK(2) 

YBAR(BHAT)  0.140515E-01  0.661724E-01  0.142004E+00  0.8850 

VARIANCE  RATIO  0 .471836E+00  0.950482E+00  0 . 999997E+00 

USING  CONTROL  VARIABLE:  WORK<3) 

YBAR(BHAT)  0.154062E-01  0.661085E-01  0.163343E+00  0.8750 

VARIANCE  RATIO  0.412782E+00  0.942611E+00  0.100000E+01 

USINC  CONTROL  VARIABLE:  W0RK(4) 

YBAR(BHAT)  0.135963E-01  0.662051E-01  0.208216E+00  0.8800 

VARIANCE  RATIO  0.471930E+00  0.947586E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR(BHAT)  0.150423E-01  0 . 7891 36E-01  0.199755E+00  0.7550 

VARIANCE  RATIO  0.282376E+00  0.849428E+00  0.999994E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.134631E-01  0.809273E-01  0.199835E+00  0.7300 

VARIANCE  RATIO  0.814486E-01  0.742618E+00  0.999995E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.1S0765E-01  0.804845E-01  0.202885E+0G  0.7450 

VARIANCE  RATIO  0.107620E+00  0.796594E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WAIT4(C/C/1) 

YBAR(BHAT)  0.148620E-01  0.810003E-01  0.200335E+00  0.7300 

VARIANCE  RATIO  0.813528E-01  0.739736E+00  0.100000E+01 
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Table  D.29.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Thirteenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 


YBAR  0 . 1 12617E+0Q  0.123006E+00  0.133489E*00 

USING  CONTROL  VARIABLE:  R0UTING(1,3) 

YBAR(BHAT)  0. 111830E+00  0 . 122966E+00  0 . 133523E+00  0.9150 

VARIANCE  RATIO  0.608003E+00  0.948123E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.112125E+00  0.123066E+00  0.135463E+00  0.8300 

VARIANCE  RATIO  0 . 395735E+00  0.839989E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0 . 1 129S8E+00  0.123095E*00  0.139967E+00  0.9050 

VARIANCE  RATIO  0.619763E+00  0.945904E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK* 2) 

YBAR(BHAT)  0.111766E+00  0.122816E+00  0.136052E+00  0.9000 

VARIANCE  RATIO  0.725878E+00  0.944308E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WORR<3) 

YBAR(BHAT)  0 . 1 13695E+00  0.123Q76E+00  0.14O833EfOO  0.9050 

VARIANCE  RATIO  0.475461E+00  0.947033E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK(4> 

YBAR(BHAT)  C.112956E+00  0.123095E+00  0.139963E+00  0.9050 

VARIANCE  RATIO  0.620128E+00  0.945910Ef00  0.999998E+00 

USING  CONTROL  VARIABLE:  SOJOURNi M/M/1 ) 

YBAR(BHAT)  0.111110E+00  0.123937E+00  0.138478Et00  0.7050 

VARIANCE  RATIO  0.165284E+00  0.690984E+00  0.998524E+00 

USING  CONTROL  VARIABLE:  WAIT4(M/M/l) 

YBAR(BHAT)  0.109029E+00  0.123774E+00  0.142480E+00  0.5100 

VARIANCE  RATIO  0.124364EfOO  0.382093E+00  0.806462E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.110603EtOO  0.123842E+00  0.13892SE+00  0.7150 

VARIANCE  RATIO  0.217436E+00  0.712158E+00  0.995266E+00 

USINC  CONTROL  VARIABLE:  WAIT4(C/G/1) 

YBAR(BHAT)  0.108567E+00  0.123744E+00  0.143333E+00  0.5050 

VARIANCE  RATIO  0 . 109572E+00  0.358926E+00  0.782032E+00 


Table  D.30.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Fourteenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERACE 


YEAR  0 . 63S066E-01  0 . 115893E<-00  0 . 172982Ef00 

USING  CONTROL  VARIABLE:  R0UTINC(1,3) 

YBAR(BHAT)  0.631326E-01  0.11566SE+00  0.174980E+00  0.9600 

VARIANCE  RATIO  0.428457E+00  0.946227E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  R0UTING<1,4) 

YBAR(BHAT)  0.634904E-01  0.116603E+00  0.207008E+00  0.9250 

VARIANCE  RATIO  0.531987E4-00  0.904762E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.639428E-01  0.11S812E+00  0.187707E+00  0.9350 

VARIANCE  RATIO  0.568546E+00  0.952948E+00  0.999991E+00 

USING  CONTROL  VARIABLE:  WORK(2) 

YBAR(BHAT)  0.627988E-01  0.116558E+00  0.206533E+00  0.9200 

VARIANCE  RATIO  0.626489E+00  0.948S41Et00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORR(3> 

YBAR(BHAT)  0.654682E-01  0.117183Ef00  0.173880E+00  0.9400 

VARIANCE  RATIO  0.629324E+00  0.956051E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK(4) 

YBAR(BHAT)  0.639314E-01  0.115819E+00  0.187862E+00  0.9350 

VARIANCE  RATIO  0.568238Ef00  0 .9529S9E*00  0.999996E+00 

USING  CONTROL  VARIABLE:  SOJOURN(M/M/l ) 

YBAR(BHAT)  0.605557E-01  0.119998E+00  0.195994E+00  0.9150 

VARIANCE  RATIO  0.371597E+00  0.922926E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR(BHAT)  0.609148E-01  0.129821E+00  0.220835E+00  0.7650 

VARIANCE  RATIO  0 . 22450SE*00  0.704063E+00  0.998094E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/C/ 1 ) 

YBAR(BHAT)  0.606823E-01  0.119896E+00  0.1959S5E+00  0.9150 

VARIANCE  RATIO  0.369715E+00  0.923695E+00  0.999985E+00 

USING  CONTROL  VARIABLE:  HAIT4(C/G/1) 

YBAR(BHAT)  0.577421E-01  0.129756E+00  0.221314E+00  0.7600 

VARIANCE  RATIO  0.226006E+00  0.700010E+00  0.997459E+00 


Table  D.31.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Fifteenth  Design  Point 


NAME  MINIMUM  AVERACE  MAXIMUM  COVERAGE 

YBAR  0 . 107831 E+00  0.132848E+00  0.155592E+00 

USING  CONTROL  VARIABLE:  ROUTINC(l,3) 

YBAR(BHAT)  0.108732E+00  0.132963E+00  0.15S574E+C0  0.9350 

VARIANCE  RATIO  0.709745E+00  0.946663E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  R0UTING(1,4) 

YBAR(BHAT)  0.106303E+00  0.132818E+00  0.155641E+00  0.9500 

VARIANCE  RATIO  0.524224Ef00  0.947768E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR(BHAT)  0.108570Ef00  0.132973E+00  0 . 1S6596E+00  0.9050 

VARIANCE  RATIO  0.751773E+00  0.949830E+00  0.999998E+00 

USING  CONTROL  VARIABLE:  WORK2) 

YBAR(BHAT)  0.8U531E-01  0.132490E+00  0.161033E+00  0.9250 

VARIANCE  RATIO  0.660564E+00  0.955374E+00  0.999989E+00 

USING  CONTROL  VARIABLE:  WORK 3) 

YBAR(BHAT)  0.999523E-01  0.132687Ef00  0.158323E+00  0.9450 

VARIANCE  RATIO  0.687302E+00  0.946703E+00  0.999967E+00 

USING  CONTROL  VARIABLE:  WORK 4) 

YBAR(BHAT)  0.108579E+00  0.132973E+00  0.156592E+00  0.9050 

VARIANCE  RATIO  0.751612E+00  0.949824E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  SOJOURN { M/M/ 1 ) 

YBAR(BHAT)  0.972244E-01  0 . 135528E*00  0.165982E+00  0.7300 

VARIANCE  RATIO  0.254447E+00  0.652769E+00  0.997814E+00 

USING  CONTROL  VARIABLE:  WAIT4( M/M/1) 

YBAR(BHAT)  0.882594E-01  0.135996E+00  0.173967E+00  0.4700 

VARIANCE  RATIO  0.586143E-01  0.2556l5Et00  0.607032E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR(BHAT)  0.884096E-01  Q.136073E+00  0.176498E+00  0.5000 

VARIANCE  RATIO  0.845022E-01  0.300370E+00  0.808305E+00 

US  INC  CONTROL  VARIABLE:  WAIT4(G/G/t) 

YBAR(BHAT)  0.895580E-01  0.135318E+00  0.173735E+00  0.4450 

VARIANCE  RATIO  0.640602E-01  0.237888E+00  0.605448E+00 


D.  32 


Table  D.32.  Control  Variate  Results  Against  Fourth  Node  Quantile 
at  the  Sixteenth  Design  Point 


NAME  MINIMUM  AVERAGE  MAXIMUM  COVERAGE 

YBAR  0.136025E-01  0.637900E-01  0.164999E+00 

USING  CONTROL  VARIABLE:  ROUTING! 1,3) 

YBAR(BHAT)  0.140895E-01  0.637S19E-01  0.163840E+00  0.8500 

VARIANCE  RATIO  0.546540E+00  0.946557E+00  0.999997E+00 

USING  CONTROL  VARIABLE:  ROUTING! 1,4) 

YBAR( BHAT)  0.133346E-01  0.632153E-01  0.165482E+00  0.8500 

VARIANCE  RATIO  0.588743E+00  0.949309E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WORK(l) 

YBAR( BHAT)  0.134627E-01  0.639593E-01  0.192337E+00  0.8550 

VARIANCE  RATIO  0.661003E+00  0.952177E+00  0.999996E+00 

USING  CONTROL  VARIABLE:  WORK! 2) 

YBAR! BHAT)  0.132382E-01  0.635254E-01  0.180854E+00  0.8500 

VARIANCE  RATIO  0.549468E+00  0.947156E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WORK! 3) 

YBAR! BHAT)  0.911578E-02  0.640720E-01  0.165083E+00  0.8300 

VARIANCE  RATIO  0.369341E+00  0.948231E+00  0.999989E+00 

USING  CONTROL  VARIABLE:  WORK! 4) 

YBAR! BHAT)  0.134616E-01  0.639375E-01  0.192181E+00  0.8550 

VARIANCE  RATIO  0.663000E+00  0 . 952478E+00  0.999970E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( M/M/ 1 ) 

YBAR! BHAT)  0.120263E-01  0.783109E-01  0.226699E+00  0.7050 

VARIANCE  RATIO  0.200864E+00  0.852099E+00  0.100000E+01 

USING  CONTROL  VARIABLE:  WAIT4(M/M/1) 

YBAR! BHAT)  -0 . 242423E-02  0.800065E-01  0.234492E+00  0.6700 

VARIANCE  RATIO  0.362815E-01  0 . 71 1689E+00  0.999992E+00 

USING  CONTROL  VARIABLE:  SOJOURN ( G/G/ 1 ) 

YBAR!  BHAT)  0 . 121 130E-01  0.8U965E-01  0 . 24561 7E+00  0.6550 

VARIANCE  RATIO  0.366392E-01  0.773461E+00  0.999999E+00 

USING  CONTROL  VARIABLE:  WAIT4!G/C/1) 

YBAR! BHAT)  -0 . 367256E-02  0.798992E-01  0.228922E+00  0.6550 

VARIANCE  RATIO  0.362530E-01  0 . 706804E+00  0.999991E+00 
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